@@ -0,0 +1,1694 @@
+%0A#!/usr/bin/python%0A%0Aimport sympy as sp%0Afrom IPython.display import display%0Aimport numpy as np%0Aimport matplotlib.pyplot as plt%0Aimport math%0Aimport matplotlib as mpl%0A%0AK = 8        %0Awg = 0.5%0Aws = 0.65%0Aalpha = 0.0%0ANA = 64%0Adef psw( j ):%0A    gc = np.log(ws/(1-ws))*(K-2*j)/(4*wg)%0A    return 0.5 + 0.5*math.erf(math.sqrt(wg)*(1.0-gc))%0A    %0A%0A# Function definition is here%0Adef tup( x ):%0A    alpha1 = alpha*(1.0-exp(-x/0.1))%0A    xx = (1.0-alpha1)*(x)%0A    return (1-x) * sum(sp.binomial(K,j) * xx**j * (1-xx)**(K-j) * psw(j) for j in xrange(0,K+1))%0A%0Adef tdown( x ):%0A    alpha1 = alpha*(1.0-exp(-x/0.1))        %0A    xx = 1.0-(1.0-alpha1)*(1.0-x)%0A    return (x) * (1 -  sum(sp.binomial(K,j) * xx**j * (1-xx)**(K-j) * psw(j) for j in xrange(0,K+1)))%0A%0Adef pswB1( x ):%0A    return  (sum(sp.binomial(K-1,j) * x**j * (1-x)**(K-1-j) * psw(j+1) for j in xrange(0,K)))%0Adef pswNB1( x ):%0A    return  (sum(sp.binomial(K-1,j) * x**j * (1-x)**(K-1-j) * psw(j) for j in xrange(0,K)))%0A    %0Adef pswB( x ):%0A    xx=alpha*(1.0+x*(float(K)-1.0))/float(K) + (1-alpha)*x%0A    return sum(sp.binomial(K,j) * xx**j * (1-xx)**(K-j) * psw(j) for j in xrange(0,K+1))%0Adef pswNB( x ):%0A    xx=alpha*(x*(float(K)-1.0))/float(K) + (1-alpha)*x%0A    return sum(sp.binomial(K,j) * xx**j * (1-xx)**(K-j) * psw(j) for j in xrange(0,K+1))%0Adef pswBNB( x ):%0A    return sum(sp.binomial(K,j) * x**j * (1-x)**(K-j) * psw(j) for j in xrange(0,K+1))%0A%0AnumX = 100%0AXs = np.zeros(numX+1)%0APAB = np.zeros(numX+1)%0Afor fX in range(0,numX+1):%0A    X = fX/float(numX)%0A    alpha = 0.8%0A    pInt = alpha%0A    f = (pInt*pswB(X) + (1.0-pInt)*pswBNB(X))/(pInt*pswNB(X) + (1.0-pInt)*pswBNB(X))%0A    Xs%5BfX%5D = X%0A    PAB%5BfX%5D = f*(X)/(f*X+(1.0-X))%0A    %0Aplt.plot(Xs,PAB)%0A%0A%0A    %0A    

test of conditional probability rate