<jupyter_start><jupyter_text>Homework 3 - Berkeley STAT 157Handout 2/5/2019, due 2/12/2019 by 4pm in Git by committing to your repository.**Formatting: please include both a .ipynb and .pdf file in your homework submission, named homework3.ipynb and homework3.pdf. You can export your notebook to a pdf either by File -> Download as -> PDF via Latex (you may need Latex installed), or by simply printing to a pdf from your browser (you may want to do File -> Print Preview in jupyter first). Please don't change the filename.**<jupyter_code>from mxnet import nd, autograd, gluon
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>1. Logistic Regression for Binary ClassificationIn multiclass classification we typically use the exponential model $$p(y|\mathbf{o}) = \mathrm{softmax}(\mathbf{o})_y = \frac{\exp(o_y)}{\sum_{y'} \exp(o_{y'})}$$1.1. Show that this parametrization has a spurious degree of freedom. That is, show that both $\mathbf{o}$ and $\mathbf{o} + c$ with $c \in \mathbb{R}$ lead to the same probability estimate.1.2. For binary classification, i.e. whenever we have only two classes $\{-1, 1\}$, we can arbitrarily set $o_{-1} = 0$. Using the shorthand $o = o_1$ show that this is equivalent to $$p(y=1|o) = \frac{1}{1 + \exp(-o)}$$1.3. Show that the log-likelihood loss (often called logistic loss) for labels $y \in \{-1, 1\}$ is thus given by $$-\log p(y|o) = \log (1 + \exp(-y \cdot o))$$1.4. Show that for $y = 1$ the logistic loss asymptotes to $o$ for $o \to \infty$ and to $\exp(o)$ for $o \to -\infty$.<jupyter_code>#1.1: C is a constant rather than the parameter in question. The model therefore does not consider c in the estiamte of the probability for the vector in question. The softmax function would end up with 1 degree of freedom due to only having 1 parameter.
#1.2: In a two class problem, a softmax function with two outputs is equivalent to a single binary output assuming a logistic function. 
# the softmax equation: p(y=1|θ)=1/1+exp(θTx);p(w2|x,θ)=1−1/1+exp(θTx)=exp(θTx)/1+exp(θTx). For the second binary function, when o-1 is set to 0, the two equations are equivalent.
#1.3: Using the linear model where y = wtx and the probability of the odds of an outcome being p/(1-p), we can
# compute the natural log as log(1/1-p). Take the natural exponential to make e^y+1 = 1/1-p. Rearranging for p gives p =1/(1+e^-y). From there, the log likelihood for p(y|o) turns into -logp(y|o) = log(1+exp(-y*o)).
#1.4 Plug in a really big number for o to test for infinity, and p(y=1|o) goes towards 1. For negative infinity, 1/(1+exp(∞)) goes towards −∞.<jupyter_output><empty_output><jupyter_text>2. Logistic Regression and Autograd1. Implement the binary logistic loss $l(y,o) = \log (1 + \exp(-y \cdot o))$ in Gluon1. Plot its values for $y \in \{-1, 1\}$ over the range of $o \in [-5, 5]$.1. Plot its derivative with respect to $o$ for $o \in [-5, 5]$ using 'autograd'.<jupyter_code>def loss(y,o):
    yhat= logistic(o)
    return -nd.nansum(y*nd.log(yhat) + (1-y)*nd.log(yhat))
def log (z):
    return 1. / (1/ +nd.exp(-z))
x = nd.arange(-5, 5, .1)
y= logistic(x)

plt.plot(x.asnumpy(),y.asnumpy())
plt.show()

def evaluate_accuracy(data_iterator, net):
    numerator = 0.
    denominator = 0.
    for i, (data, label) in enumerate(data_iterator):
        data = data.as_in_context(model_ctx).reshape((-1,784))
        label = label.as_in_context(model_ctx)
        label_one_hot = nd.one_hot(label, 10)
        output = net(data)
        predictions = nd.argmax(output, axis=1)
        numerator += nd.sum(predictions == label)
        denominator += data.shape[0]
    return (numerator / denominator).asscalar()
epochs = 10
loss_sequence = []
num_examples = len(Xtrain)

for e in range(epochs):
    cumulative_loss = 0
    for i, (data, label) in enumerate(train_data):
        data = data.as_in_context(model_ctx)
        label = label.as_in_context(model_ctx)
        with autograd.record():
            output = net(data)
            loss = log_loss(o, label)
        loss.backward()
        trainer.step(batch_size)
        cumulative_loss += nd.sum(loss).asscalar()
    print("Epoch %s, loss: %s" % (e, cumulative_loss ))
    loss_sequence.append(cumulative_loss)
    
%matplotlib inline

import matplotlib
import matplotlib.pyplot as plt

plt.figure(num=None,figsize=(10, 10))
plt.plot(loss_sequence)<jupyter_output><empty_output><jupyter_text>3. Ohm's LawImagine that you're a young physicist, maybe named [<NAME>](https://en.wikipedia.org/wiki/Georg_Ohm), trying to figure out how current and voltage depend on each other for resistors. You have some idea but you aren't quite sure yet whether the dependence is linear or quadratic. So you take some measurements, conveniently given to you as 'ndarrays' in Python. They are indicated by 'current' and 'voltage'.Your goal is to use least mean squares regression to identify the coefficients for the following three models using automatic differentiation and least mean squares regression. The three models are:1. Quadratic model where $\mathrm{voltage} = c + r \cdot \mathrm{current} + q \cdot \mathrm{current}^2$.1. Linear model where $\mathrm{voltage} = c + r \cdot \mathrm{current}$.1. Ohm's law where $\mathrm{voltage} = r \cdot \mathrm{current}$.<jupyter_code>current = nd.array([1.5420291, 1.8935232, 2.1603365, 2.5381863, 2.893443, \
                    3.838855, 3.925425, 4.2233696, 4.235571, 4.273397, \
                    4.9332876, 6.4704757, 6.517571, 6.87826, 7.0009003, \
                    7.035741, 7.278681, 7.7561755, 9.121138, 9.728281])
voltage = nd.array([63.802246, 80.036026, 91.4903, 108.28776, 122.781975, \
                    161.36314, 166.50816, 176.16772, 180.29395, 179.09758, \
                    206.21027, 272.71857, 272.24033, 289.54745, 293.8488, \
                    295.2281, 306.62274, 327.93243, 383.16296, 408.65967])
linear_model = linreg(current, voltage, c)
Ohms = linreg(current, volt)


w = nd.random_normal(shape=(1,1), ctx=model_ctx)
r = nd.random_normal(shape=(1,1), ctx=model_ctx)
c = nd.random_normal(shape = (1,1), ctx = model_ctx)
q = nd.random_normal(shape = (1,1), ctx = model_ctx)

params = [w, r, c, q]
for param in params:
    param.attach_grad()
def Ohms(x):
    return x * r
def voltage(r):
    return mx.nd.dot(r, w)+ c
def quad_voltage(X):
    return r * w + q * w^2 + c<jupyter_output><empty_output><jupyter_text>4. EntropyLet's compute the *binary* entropy of a number of interesting data sources. 1. Assume that you're watching the output generated by a [monkey at a typewriter](https://en.wikipedia.org/wiki/File:Chimpanzee_seated_at_typewriter.jpg). The monkey presses any of the $44$ keys of the typewriter at random (you can assume that it has not discovered any special keys or the shift key yet). How many bits of randomness per character do you observe?1. Unhappy with the monkey you replaced it by a drunk typesetter. It is able to generate words, albeit not coherently. Instead, it picks a random word out of a vocabulary of $2,000$ words. Moreover, assume that the average length of a word is $4.5$ letters in English. How many bits of randomness do you observe now?1. Still unhappy with the result you replace the typesetter by a high quality language model. These can obtain perplexity numbers as low as 20 points per character. The perplexity is defined as a length normalized probability, i.e.$$\mathrm{PPL}(x) = \left[p(x)\right]^{1/\mathrm{length}(x)}$$<jupyter_code>#1. 6 bits of randomness per character
#2. 11.8 bits per word, 2.6 bits per letter<jupyter_output><empty_output><jupyter_text>5. Wien's Approximation for the Temperature (bonus)We will now abuse Gluon to estimate the temperature of a black body. The energy emanated from a black body is given by Wien's approximation.$$B_\lambda(T) = \frac{2 h c^2}{\lambda^5} \exp\left(-\frac{h c}{\lambda k T}\right)$$That is, the amount of energy depends on the fifth power of the wavelength $\lambda$ and the temperature $T$ of the body. The latter ensures a cutoff beyond a temperature-characteristic peak. Let us define this and plot it.<jupyter_code># Lightspeed
c = 299792458
# Planck's constant
h = 6.62607004e-34
# Boltzmann constant
k = 1.38064852e-23
# Wavelength scale (nanometers)
lamscale = 1e-6
# Pulling out all powers of 10 upfront
p_out = 2 * h * c**2 / lamscale**5
p_in = (h / k) * (c/lamscale) 

# Wien's law
def wien(lam, t):
    return (p_out / lam**5) * nd.exp(-p_in / (lam * t))

# Plot the radiance for a few different temperatures
lam = nd.arange(0,100,0.01)
for t in [10, 100, 150, 200, 250, 300, 350]:
    radiance = wien(lam, t)
    plt.plot(lam.asnumpy(), radiance.asnumpy(), label=('T=' + str(t) + 'K'))
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Next we assume that we are a fearless physicist measuring some data. Of course, we need to pretend that we don't really know the temperature. But we measure the radiation at a few wavelengths.<jupyter_code># real temperature is approximately 0C
realtemp = 273
# we observe at 3000nm up to 20,000nm wavelength
wavelengths = nd.arange(3,20,2)
# our infrared filters are pretty lousy ...
delta = nd.random_normal(shape=(len(wavelengths))) * 1

radiance = wien(wavelengths + delta,realtemp)
plt.plot(wavelengths.asnumpy(), radiance.asnumpy(), label='measured')
plt.plot(wavelengths.asnumpy(), wien(wavelengths, realtemp).asnumpy(), label='true')
plt.legend()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12313551253570296, 'fraction_numerical': 0.05691314926478366, 'mean_word_length': 3.9575249082328265, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 33, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '308360', 'license': 'Apache-2.0', 'n_tokens_mistral': 3529, 'n_tokens_neox': 3104, 'n_words': 1132, 'path': 'homeworks/homework3 (2).ipynb', 'repo_name': 'SWChhh20/berkeley-stat-157', 'size': '9462'}