#!/usr/bin/env python
# -*- coding: utf-8 -*-

# ***************************************************************************
# *   Copyright (C) 2014 by Paul Lutus                                      *
# *   lutusp@arachnoid.com                                                  *
# *                                                                         *
# *   This program 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 2 of the License, or     *
# *   (at your option) any later version.                                   *
# *                                                                         *
# *   This program 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 program; if not, write to the                         *
# *   Free Software Foundation, Inc.,                                       *
# *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
# ***************************************************************************

# Created on Feb 26, 2014 8:47:39 AM

import re,sys,os,math

def f(x,m,sd):
  sq2 = math.sqrt(2.0)
  # see Sage worksheet 'p-value_sigma_erf_cdf'
  return math.erf((sq2 * (x - m))/(2 * sd)) * .5
  
def int_f(a,b,m,sd):
  return f(b,m,sd) - f(a,m,sd)
  
def draw_line(n):
  print('%s' % '-' * n)

if(len(sys.argv) == 5):
  (a,b,m,sd) = (float(x) for x in sys.argv[1:])
  print('Equation: y = ½ (erf((√2 (b-μ)) / 2 σ) - erf((√2 (a-μ)) / 2 σ))')
  print('Arguments:')
  print('%12s %12s %13s %13s' % ('a','b','μ','σ'))
  draw_line(51)
  print('%12.8g %12.8g %12.8g %12.8g' % (a,b,m,sd))
  print('Result:')
  print('%12s %12s' % ('y','1-y'))
  draw_line(25)
  y = int_f(a,b,m,sd)
  print('%12.8g %12.8g' % (y,1.0-y))
else:
  print('usage: %s lowerbound upperbound μ σ' % sys.argv[0])
  print('%13s %22s %22s %22s %22s' % ('σ','two-tailed: y','two-tailed: 1-y','one-tailed: y','one-tailed: 1-y'))
  for x in range(1,10):
    a = int_f(-x,x,0,1)
    b = int_f(-1e9,x,0,1)
    print('%12f %22.16g %22.16g %22.16g %22.16g' % (x,a,1.0-a,b,1.0-b))