<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Precision vs Accuracy Precision is the accuracy of basic arithmetic operations used in the computation Accuracy is the absolute or relative error of the approximate quantity* NOTE
Step1:  The last definition avoids issues that the previous two had, but it's still non-ideal
Step2:  Absolute and relative error
Step3:  Relative Error for Non-ScalarsFor non-scalars, calculating normalized value ||x|| = max / sum, implies smaller components of x are bound by absolute error only. Consider compensative relative error
Step4:  Measuring Accuracy## Backward and Forward Errors Forward error - error of the computed value Backward error
Step5:  Quoting Wikipedia
Step6:  Rounding Multiple Times Can Accumulate Error
Step7:  (Optional) Differences between round and np.roundThere are few differences between built-in Python 2.7 round function and numpy (a)round
Step8:  Loss of significanceError in floating point arithmetic when an operation increases relative error substantially more than absolute error
Step9:  CancellationCancellation is an example of loss of significance and it happens when two nearly equal numbers are subtracted and can lead to significant inaccuracies. As an example let's look at function
Step10:  Fortunately, we can rewrite f(x) in a form that is less prone to cancellation
Step11:  Summing numbers
Step12:  Calculations without subtractions are fine, right?Hint
Step13:  Sometimes the rounding errors cancel out and produce a result more accurate than the intermediate calculations
Step14:  Rounding errors are not random
Step15:  Examples## Variance CalculationThere are two concepts that we refer to as variance
<ASSISTANT_TASK:>
Python Code:
# Definition 1: Round down to p-sig. digit numberx  = 0.90x1 = 0.99 # 2 correct significant digit, actual difference 0.09x2 = 0.89 # 1 correct significant digit, actual difference 0.01
# Definition 2: Round to nearest p-sig. digit numbery  = 0.9951 # --> 0.10y1 = 0.9499 # --> 0.90 , only 1 correct sig. digity2 = 1.0000 # --> 0.10 , 3 correct sig. digits
Explanation: Precision vs Accuracy Precision is the accuracy of basic arithmetic operations used in the computation Accuracy is the absolute or relative error of the approximate quantity* NOTE: Accuracy is not limited by precision, finite precision arithmetic can simulate any precision with more computation# Measuring Precision## Significant DigitsThe number of significant digits may be imprecise, prefer to use relative error.
End of explanation
# Definition 3: Numbers x and x' match to p-sig. digits if x - x' < half a unit in p-th sig. digit of xx1 = 0.123x2 = 0.127# 0.004 < (0.01 / 2) => x1 and x2 match in 2 significant digits according to this definition wchich may be slightly confusing 
Explanation: The last definition avoids issues that the previous two had, but it's still non-ideal
End of explanation
def absolute_error(true_value, approx_value):    return abs(true_value - approx_value)print 'Absolute error: {0:.9f}'.format(absolute_error(10.951, 10.949))def relative_error(true_value, approx_value):    return absolute_error(true_value, approx_value) / abs(true_value)print 'Relative error: {0:.9f}'.format(relative_error(10.951, 10.949))
Explanation: Absolute and relative error
End of explanation
import numpy as npdef relative_error(true_value, approx_value):    return np.max(np.fabs((true_value - approx_value) / true_value))x_value = np.array([10000, 0.01])x_approx = np.array([9999, 0.02])                   print relative_error(x_value, x_approx)
Explanation: Relative Error for Non-ScalarsFor non-scalars, calculating normalized value ||x|| = max / sum, implies smaller components of x are bound by absolute error only. Consider compensative relative error: max(i) |xi - xi’| / |xi|, which puts all components on equal footing.
End of explanation
# Wilkinson's polynomial is defined as p(x) = (x - 1)(x - 2)...(x - 20)import numpy as npx = np.linspace(0, 20, 4000)y = 1for i in range(1, 20):    y *= (x - i)import matplotlib.pyplot as plt%matplotlib inlineplt.ylim([-4e11, 4e11])plt.plot(x, y)
Explanation: Measuring Accuracy## Backward and Forward Errors Forward error - error of the computed value Backward error:  * Let y = f(x), given x we approximate f(x) with y’  * Let dx be the smallest quantity where y’ = f(x + dx) in exact computation  * Then dx is the backward error Benefits of using backward error:  * Unifies error w/ perturbation in the data  * Removes the need to calculate forward error Forward-backward error:  * f(x + dx) = y + dy  * Used to define stability of computation where just using backward error isn’t possible, e.g. sin, cos If rounding errors are dominant source of errors, we call an algorithm numerically stable if it is stable in forward-backward error sense## Condition Number Condition number of a function with respect to its arguments is used to measure how how much the output of the function will change for a small change in the input As a rule of thumb, if the condition number kappa(A) = 10^k, then you may lose up to k digits of accuracy on top of what would be lost to the numerical method due to loss of precision from arithmetic methods A problem with a low condition number is said to be well-conditioned, while a problem with a high condition number is said to be ill-conditioned For example, the condition number associated with the linear equation Ax = b gives a bound on how inaccurate the solution x will be after approximation. Note that this is before the effects of round-off error are taken into account; conditioning is a property of the matrix, not the algorithm or floating point accuracy of the computer used to solve the corresponding system.An important example of an ill-conditioned problem is finding roots of a polynomial*. Let's look at Wilkinson’s polynomial as an example.
End of explanation
import numpy as npx = np.float64(0.1)y = np.float32(0.1)print '32- vs 64-bit representation difference:', abs(x - y)
Explanation: Quoting Wikipedia: "If the coefficient of x19 is decreased from −210 by 2−23 to −210.0000001192, then the polynomial value w(20) decreases from 0 to −6.25×10^17, and the root at x = 20 grows to x ≈ 20.8. The roots at x = 18 and x = 19 collide into a double root at x ≈ 18.62 which turns into a pair of complex conjugate roots at x ≈ 19.5±1.9i as the perturbation increases further." *"Wilkinson's polynomial is often used to illustrate the undesirability of naively computing eigenvalues of a matrix by first calculating the coefficients of the matrix's characteristic polynomial* and then finding its roots, since using the coefficients as an intermediate step may introduce an extreme ill-conditioning even if the original problem was well conditioned"
Sources of Error## Truncation error (discretization error)Error coming from representing a function or continuous variable using finite number of evaluations - outside of scope of this notebook, mentioned for completeness only.## Round-off error Difference between calculated approximation and exact value due to rounding Related to representation error, which is due to representing numbers with finite number of digits
End of explanation
import numpy as np# For explanation why we we use np.round rather than default Python 2.7.3 round function, see belowx = 9.945309 print np.round(x, 2), np.round(np.round(x, 2), 1)print np.round(x, 1)
Explanation: Rounding Multiple Times Can Accumulate Error
End of explanation
import numpy as npfor i in range(13):    x = -3 + 0.5 * i    print '\t{0:5.1f}\t{1:5.1f}\t{2:5.1f}'.format(x, round(x), np.round(x))
Explanation: (Optional) Differences between round and np.roundThere are few differences between built-in Python 2.7 round function and numpy (a)round: The built in function rounds away from zero Numpy round rounds to even, which tends to skew the results less and is a commonly accepted rounding method* From my (limited) exerience it looks numpy round is much better behaved in dealing with decimal-to-binary float rounding errorsNote that Python 3 has a different round function that behaves more similarly to numpy round.
End of explanation
x_value = 0.123123y_value = 0.123000# We want to learn the value of d = x - y:d_value = x_value - y_valueprint 'Actual d value:', # Assuming we're apprximating above calculation with a 4 decimal digits precision:import numpy as npd_approx = np.round(x_value, 4) - np.round(y_value, 4)print 'Approx d value:', d_approxprint 'Absolute error: {0:.9f}'.format(abs(d_value - d_approx))print 'Relative error: {0:.9f}'.format(abs((d_value - d_approx) / d_value))
Explanation: Loss of significanceError in floating point arithmetic when an operation increases relative error substantially more than absolute error
End of explanation
import numpy as npfrom math import sin, cosnear_zero = 1.2e-8def f(x):    return (1 - np.cos(x)) / x**2print 'Value of f near zero:', f(near_zero)x = np.linspace(near_zero, 10, 100)y = f(x)import matplotlib.pyplot as plt%matplotlib inlineplt.ylim([0, 1])plt.plot(x, y)
Explanation: CancellationCancellation is an example of loss of significance and it happens when two nearly equal numbers are subtracted and can lead to significant inaccuracies. As an example let's look at function:  f(x) = (1 - cos x) / x^2  We can see that calculating f(x) near zero may lead to issues, since cos(0) = 1.
End of explanation
def g(x):    return 0.5 * (2 * np.sin(x / 2) / x)**2print 'Value of g near zero:', g(near_zero)x = np.linspace(near_zero, 10, 100)y = g(x)import matplotlib.pyplot as plt%matplotlib inlineplt.ylim([0, 1])plt.plot(x, y)
Explanation: Fortunately, we can rewrite f(x) in a form that is less prone to cancellation:
End of explanation
from math import fsumprint '{0:0.20f}'.format(sum([0.1] * 10))print '{0:0.20f}'.format(fsum([0.1] * 10))
Explanation: Summing numbers: https://docs.python.org/2/library/math.html https://en.m.wikipedia.org/wiki/Kahan_summation_algorithmInterestingly, IPython Notebook appears to do the right thing by default, what?IPython session:In [6]: from math import fsumIn [7]: sum([0.1] * 10) Out[7]: 0.9999999999999999In [8]: fsum([0.1] * 10) Out[8]: 1.0
End of explanation
def naive_e(n):    return (1 + 1.0 / n)**nfrom math import expe = exp(1)for i in range(5, 20):    print naive_e(10**i) - e
from math import sqrtdef identity(x, n):    for i in xrange(n):        x = sqrt(x)    for i in xrange(n):        x = x**2    return xx = 2for i in xrange(35, 60):    print x - identity(x, i)
Explanation: Calculations without subtractions are fine, right?Hint: No, we can still over/under-flow the floating point precision numbers.
End of explanation
# Computing f(x) = (exp(x) - 1) / x == sum(x^i / (i + 1)!)from math import exp, logdef f1(x):    if 0 == x:        return 1    return (exp(x) - 1) / xdef f2(x):    if 0 == x:        return 1    y = exp(x)    return (y - 1) / log(y)# f(epsilon) ~= 1for i in range(8, 15):    epsilon = 1.0 / (10**i)    print 'epsilon:', epsilon    print '|1 - f1(epsilon)|:', abs(1 - f1(epsilon))    print '|1 - f2(epsilon)|:', abs(1 - f2(epsilon))    print# NOTE: Above doesn't hold if we calculate for powers of 2!# for i in range(30, 40):#     print abs(1 - f1(1.0 / (2**i))), abs(1 - f2(1.0 / (2**i)))
Explanation: Sometimes the rounding errors cancel out and produce a result more accurate than the intermediate calculations
End of explanation
import numpy as npdef r(x):    # Calculate value of a rational function using Horner's rule (https://en.wikipedia.org/wiki/Horner%27s_method)    # More on the function can be looked up on Wolfram's Alpha:    # http://www.wolframalpha.com/input/?i=f(x)+%3D+(622.0+-+x+*+(751.0+-+x+*+(324.0+-+x+*+(59.0+-+4+*+x))))+%2F+(112+-+x+*+(151+-+x+*+(72+-+x+*+(14+-+x))))&t=crmtb01    p = 622.0 - x * (751.0 - x * (324.0 - x * (59.0 - 4 * x)))    q = 112 - x * (151 - x * (72 - x * (14 - x)))    return p / qdef calc_f(a):    t = np.array([a + k * 2**-52 for k in xrange(400)])    t = r(t)    t -= t[0]    t *= 1.0 / max(abs(t))    return timport matplotlib.pyplot as plt%matplotlib inlinedef plot(t):    plt.plot(t, linestyle='--', marker='o')    plt.show()    for a in [1.606, 4, 8, 16, 32]:    plot(calc_f(a))
Explanation: Rounding errors are not random
End of explanation
# Calculate the variance of X = sum of dice in N throws. The probabilities of each side of dice are given.import numpy as np# Dice sides probabilities, p[i] = probability of throwing i. We use fractions to be able to compare our calculations to# 'exact' values:from fractions import Fraction as Ffrom numpy.random import randint as rimax = 1000000pf = np.array([0] + [F(ri(max), 1) for _ in range(6)])pf /= sum(pf)p = np.array([float(f) for f in pf])# Number of throwsN = 10000 # 30000# Dynamic program, we're holding probability of throwing k-dice for each possibility.# First iteration:dp = np.ones(1)for i in range(N):    dp = np.convolve(dp, p)print dp# Let's calculate variance using both ways.dice = np.arange(len(dp))# var1 = E[(X - E[X])^2]ex = (p * np.arange(7)).sum() * Nvar1 = (dp * (dice - ex)**2).sum()print 'Variance calculated using: E[(X - E[X])^2]:   ', var1# var2 = E[X^2] - (E[X])^2ex2 = (dice**2 * dp).sum()var2 = ex2 - ex**2print 'Variance calculated using: E[X^2] - (E[X])^2: ', var2print 'E[X^2]:    ', ex2print '(E[X2])^2: ', ex**2# There is a simpler way to calculate variance in this particular problem, since all N throws are independent we can simply# calculate variability of one and multiply it by N. To make sure this calculation is as precise as possible, we use Python# representation of rational numbers.# Let's use equation 1: E[(X - E[X])^2]rangef = np.array([F(i, 1) for i in range(7)])exf = (rangef * pf).sum()ex2f = (rangef**2 * pf).sum()varf = ex2f - exf**2varf *= Nprint 'Variance calculated fractions:                ', float(varf)print 'Relative error of: E[(X - E[X])^2]:   ', abs((var1 - varf) / varf)print 'Relative error of: E[X^2] - (E[X])^2: ', abs((var2 - varf) / varf)
Explanation: Examples## Variance CalculationThere are two concepts that we refer to as variance:0. A property of a distribution0. A characteristic of a set of observationsVariance of a random variable X is defined as:  Var(X) = E[(X - E[X])^2] = ... = E[X^2] - (E[X])^2  Where E[Z] is the expected value of a random variable Z.The second form (E[X^2] - (E[X])^2) should be avoided when performing calculations on a fixed precision machine. Although it has a nice property that for sample variance it is easy to implement it while traversing the data just once, naive implementations usually suffer from extreme cancellation.
End of explanation
<END_TASK>