#!/usr/bin/env python
#-*- coding:utf-8 -*-
"""
This file is part of www.cogsci.nl.
For more info, see
This is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This software is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this. If not, see .
"""
from math import *
import random
import bayes
import numpy
import scipy.stats
def dnorm(x, M = 0, var = 1):
"""
Returns the height of the normal distribution at a given point
Arguments:
x -- a point along the x-axis
Keyword arguments:
M -- the mean of the distribution
var -- the variance (sd ** 2) of the distribution
Returns:
The height of the normal distribution at point x
"""
return (1.0 / sqrt(2.0 * pi * var)) * e ** -( (x - M) ** 2.0 / (2.0 * var))
def uniform_bf(M, SE, lower = 0, upper = 1):
"""
Calculates a uniform Bayes factor that estimates the evidence over
the null hypothesis. Based on the R-code by Baguley & Kaye (2010).
Arguments:
M -- the observed mean
SE -- the standard error of the measurement
Keyword arguments:
lower -- a lower bound of the predicted value (default = 0)
upper - an upper bound of the predicted value (upper = 1)
Returns:
A dictionary containing 'likelyhood_theory', 'likelyhood_null' and 'bf' (Bayes factor)
"""
area = 0
theta = lower
_range = upper - lower
incr = _range / 2000
for A in range(-1000, 1000):
theta += incr
dist_theta = 1 / _range
height = dist_theta * dnorm(M, M = theta, var = SE ** 2)
area += height * incr
likelyhood_theory = area
likelyhood_null = dnorm(M, M = 0, var = SE ** 2)
bf = likelyhood_theory / likelyhood_null
return {"likelyhood_theory" : likelyhood_theory, "likelyhood_null" : likelyhood_null, "bf" : bf}
random.seed()
# A dictionary to keep track of false alarms
bf_fa = {}
p_fa = {}
runs = 500 # Nr of simulated experiments
samples = 5000 # Nr of subjects
# Walk through all 'experiments'
for run in range(runs):
data = [] # A list of datapoints
p_hit = False # Has a p-value false alarm occurred?
bf_hit = False # Has a bayesian false alarm occurred?
# Walk through all 'subjects'
for i in range(samples):
# Get a random datapoint and store it in the datalist
val = random.random() - 0.5
data.append(val)
# If we have enough data to start testing (at least two)
if len(data) > 2:
# Calculate the bayes factor if a false alarm hasn't occurred yet (afterwards it doesn't matter)
if not bf_hit:
bf = uniform_bf(numpy.mean(data), numpy.std(data) / sqrt(len(data)), lower = -0.5, upper = 0.5)["bf"]
# Calculate the p-value if a false alarm hasn't occurred yet (afterwards it doesn't matter)
if not p_hit:
t, p = scipy.stats.ttest_1samp(data, 0)
# Make sure that the dictionary contains an entry for this subject nr
if i not in p_fa:
p_fa[i] = 0
if i not in bf_fa:
bf_fa[i] = 0
# Determine a p-value false alarm and remember it
if p < .05:
p_hit = True
if p_hit:
p_fa[i] += 1
# Determine a bayesian false alarm and remember it
if bf > 3:
bf_hit = True
if bf_hit:
bf_fa[i] += 1
# Output the results
for i in range(2, samples):
print "%d\t%d\t%d" % (i, bf_fa[i], p_fa[i])