# Scalars


```python
%matplotlib inline
```


```python
import numpy as np
import matplotlib.pyplot as plt
```

## Integers

### Binary representation of integers


```python
format(16, '032b')
```




    '00000000000000000000000000010000'



### Bit shifting


```python
format(16 >> 2, '032b')
```




    '00000000000000000000000000000100'




```python
16 >> 2
```




    4




```python
format(16 << 2, '032b')
```




    '00000000000000000000000001000000'




```python
16 << 2
```




    64



### Overflow

In general, the computer representation of integers has a limited range, and may overflow. The range depends on whether the integer is signed or unsigned.

For example, with 8 bits, we can represent at most $2^8 = 256$ integers.

- 0 to 255 unsigned
- -128 ti 127 signed

Signed integers


```python
np.arange(130, dtype=np.int8)[-5:]
```




    array([ 125,  126,  127, -128, -127], dtype=int8)



Unsigned integers


```python
np.arange(130, dtype=np.uint8)[-5:]
```




    array([125, 126, 127, 128, 129], dtype=uint8)




```python
np.arange(260, dtype=np.uint8)[-5:]
```




    array([255,   0,   1,   2,   3], dtype=uint8)



### Integer division

In Python 2 or other languages such as C/C++, be very careful when dividing as the division operator `/` performs integer division when both numerator and denominator are integers. This is rarely what you want. In Python 3 the `/` always performs floating point division, and you use `//` for integer division, removing a common source of bugs in numerical calculations.


```python
%%python2

import numpy as np

x = np.arange(10)
print(x/10)
```

    [0 0 0 0 0 0 0 0 0 0]


Python 3 does the "right" thing.


```python
x = np.arange(10)
x/10
```




    array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])



## Real numbers

Real numbers are represented as **floating point** numbers. A floating point number is stored in 3 pieces (sign bit, exponent, mantissa) so that every float is represetned as get +/- mantissa ^ exponent. Because of this, the interval between consecutive numbers is smallest (high precison) for numebrs close to 0 and largest for numbers close to the lower and upper bounds.

Because exponents have to be singed to represent both small and large numbers, but it is more convenint to use unsigned numbers here, the exponnent has an offset (also knwnn as the exponentn bias). For example, if the expoennt is an unsigned 8-bit number, it can rerpesent the range (0, 255). By using an offset of 128, it will now represent the range (-127, 128).



**Note**: Intervals between consecutive floating point numbers are not constant. In particular, the precision for small numbers is much larger than for large numbers. In fact, approximately half of all floating point numbers lie between -1 and 1 when using the `double` type in C/C++ (also the default for `numpy`).



Because of this, if you are adding many numbers, it is more accurate to first add the small numbers before the large numbers.

#### IEEE 754 32-bit floating point representation



See [Wikipedia](https://en.wikipedia.org/wiki/Single-precision_floating-point_format) for how this binary number is evaluated to 0.15625.


```python
from ctypes import c_int, c_float
```


```python
s = c_int.from_buffer(c_float(0.15625)).value
```


```python
s = format(s, '032b')
s
```




    '00111110001000000000000000000000'




```python
rep = {
    'sign': s[:1], 
    'exponent' : s[1:9:], 
    'fraction' : s[9:]
}
rep
```




    {'sign': '0', 'exponent': '01111100', 'fraction': '01000000000000000000000'}



### Most base 10 real numbers are approximations

This is simply because numbers are stored in finite-precision binary format.


```python
'%.20f' % (0.1 * 0.1 * 100)
```




    '1.00000000000000022204'



### Never check for equality of floating point numbers


```python
i = 0
loops = 0
while i != 1:
    i += 0.1 * 0.1
    loops += 1
    if loops == 1000000:
        break
i
```




    10000.000000171856




```python
i = 0
loops = 0
while np.abs(1 - i) > 1e-6:
    i += 0.1 * 0.1
    loops += 1
    if loops == 1000000:
        break
i
```




    1.0000000000000007



### Associative law does not necessarily hold


```python
6.022e23 - 6.022e23 + 1
```




    1.0




```python
1 + 6.022e23 - 6.022e23
```




    0.0



### Distributive law does not hold


```python
a = np.exp(1)
b = np.pi
c = np.sin(1)
```


```python
a*(b+c)
```




    10.82708950985241




```python
a*b + a*c
```




    10.827089509852408



### Catastrophic cancellation

Consider calculating sample variance

$$
s^2= \frac{1}{n(n-1)}\sum_{i=1}^n x_i^2 - (\sum_{i=1}^n x_i)^2
$$

Be  careful whenever you calculate the difference of potentially big numbers.


```python
def var(x):
    """Returns variance of sample data using sum of squares formula."""
    
    n = len(x)
    return (1.0/(n*(n-1))*(n*np.sum(x**2) - (np.sum(x))**2))
```

### Underflow


```python
np.warnings.filterwarnings('ignore')
```


```python
np.random.seed(4)
xs = np.random.random(1000)
ys = np.random.random(1000)
np.prod(xs)/np.prod(ys)
```




    nan



#### Prevent underflow by staying in log space


```python
x = np.sum(np.log(xs))
y = np.sum(np.log(ys))
np.exp(x - y)
```




    696432868222.2549



### Overflow


```python
np.exp(1000)
```




    inf



### Numerically stable algorithms

#### What is the sample variance for numbers from a normal distribution with variance 1?


```python
np.random.seed(15)
x_ = np.random.normal(0, 1, int(1e6))
x = 1e12 + x_
var(x)
```




    -1328166901.4739892



#### Use functions from numerical libraries where available


```python
np.var(x)
```




    1.001345504504934



There is also a variance function in the standard library, but it is slower for large arrays.


```python
import statistics

statistics.variance(x)
```




    1.0013464175061981



Note that `numpy` uses does not use the asymptotically unbiased estimator by default. If you want the unbiased variance, set `ddof` to 1.


```python
np.var([1,2,3,4], ddof=1)
```




    1.6666666666666667




```python
statistics.variance([1,2,3,4])
```




    1.6666666666666667



### Useful numerically stable functions 

Let's calculate

$$
\log(e^{1000} + e^{1000})
$$

Using basic algebra, we get the solution $\log(2) + 1000$.

\begin{align}
\log(e^{1000} + e^{1000}) &= \log(e^{0}e^{1000} + e^{0}e^{1000}) \\
&= \log(e^{100}(e^{0} + e^{0})) \\
&= \log(e^{1000}) + \log(e^{0} + e^{0}) \\
&= 1000 + \log(2)
\end{align}

**logaddexp**


```python
x = np.array([1000, 1000])
np.log(np.sum(np.exp(x)))
```




    inf




```python
np.logaddexp(*x)
```




    1000.6931471805599



**logsumexp**

This function generalizes `logaddexp` to an arbitrary number of addends and is useful in a variety of statistical contexts.

Suppose we need to calculate a probability distribution $\pi$ parameterized by a vector $x$

$$
\pi_i = \frac{e^{x_i}}{\sum_{j=1}^n e^{x_j}}
$$

Taking logs, we get

$$
\log(\pi_i) = x_i - \log{\sum_{j=1}^n e^{x_j}}
$$


```python
x = 1e6*np.random.random(100)
```


```python
np.log(np.sum(np.exp(x))) 
```




    inf




```python
from scipy.special import logsumexp
```


```python
logsumexp(x)
```




    984755.5793729819



**logp1 and expm1**


```python
np.exp(np.log(1 + 1e-6)) - 1
```




    9.999999999177334e-07




```python
np.expm1(np.log1p(1e-6))
```




    1e-06



**sinc**


```python
x = 1
```


```python
np.sin(x)/x
```




    0.8414709848078965




```python
np.sinc(x)
```




    3.8981718325193755e-17




```python
x = np.linspace(0.01, 2*np.pi, 100)
```


```python
plt.plot(x, np.sinc(x), label='Library function')
plt.plot(x, np.sin(x)/x, label='DIY function')
plt.legend()
pass
```

cb54ee6e7e06d88b96e4f63564c0a2d4dcc74690
ipynb
Jupyter Notebook
notebook/S03A_Scalars_Annotated.ipynb
ashnair1/sta-663-2019
17eb85b644c52978c2ef3a53a80b7fb031360e3d
BSD-3-Clause
2019-01-09T21:53:55.000Z
2022-02-16T17:14:22.000Z
notebook/S03A_Scalars_Annotated.ipynb
ashnair1/sta-663-2019
17eb85b644c52978c2ef3a53a80b7fb031360e3d
BSD-3-Clause
notebook/S03A_Scalars_Annotated.ipynb
ashnair1/sta-663-2019
17eb85b644c52978c2ef3a53a80b7fb031360e3d
BSD-3-Clause
2019-01-09T21:43:48.000Z
2021-11-15T04:26:25.000Z
Qwen/Qwen-72B
1. YES
2. YES