Commit 92a660b9 by Erik Strand

Start pset 11

parent 6a7bb75e
 import numpy as np import math N = 100 x = np.random.uniform(-5.0, 5.0, N) y = np.tanh(x) + np.random.normal(0.0, 0.25, N) # number of clusters M = 1 # prevents variances from shrinking to zero tiny = 0.0 # initialize input distribution means and variances input_means = np.random.uniform(-5.0, 5.0, M) print("input_means") print(input_means) input_vars = np.array([ 1.0 ] * M) print("input_vars") print(input_vars) # define input distributions p(x | c_m) def input_dist(x, m): mu = input_means[m] var = input_vars[m] scale = 1 / (np.sqrt(2.0 * math.pi * var)) return scale * np.exp(-(x - mu)**2.0 / (2 * var)) print("input probs at 0") print([input_dist(0.0, i) for i in range(M)]) # initialize function parameters beta = [ np.array([0, 1]) ] * M print("beta") print(beta) # initialize functions def cluster_func(x, m): return beta[m][0] + beta[m][1] * x # initialize output variances out_var = [ 1.0 ] * M print("output variances") print(out_var) # define output distributions p(y | x, c_m) def output_dist(y, x, m): mu = cluster_func(x, m) var = out_var[m] scale = 1 / (np.sqrt(2.0 * math.pi * var)) return scale * np.exp(-(y - mu)**2.0 / (2 * var)) print("input probs for 0 -> 0") print([output_dist(0.0, 0.0, i) for i in range(M)]) # initialize cluster probabilities p(c_m) p_c = [ 1.0 / M ] * M print("cluster probs p(c_m)") print(p_c) # define posterior distributions p(c_m | y, x) # TODO denominator is shared, can calculate once each cycle def posterior(m, y, x): numerator = output_dist(y, x, m) * input_dist(x, m) * p_c[m] denominator = sum( output_dist(y, x, i) * input_dist(x, i) * p_c[i] for i in range(M)) return numerator / denominator print("posterior(cluster 0 | y = 0, x = 0)") print(posterior(0, 0.0, 0.0)) # define cluster weighted expectation helper # values should be an N long array (your quantity of interest evaluated for all samples) def cw_expectation(values, m): numerator = sum(values[n] * posterior(m, y[n], x[n]) for n in range(N)) denominator = N * p_c[m] # equivalently, could use #denominator = sum(posterior(m, y[n], x[n]) for n in range(N)) return numerator / denominator def em_step(): print("================================================================================") # update cluster probabilities # TODO posterior should operate on vectors directly for m in range(M): p_c[m] = (1.0 / N) * sum(posterior(m, y[n], x[n]) for n in range(N)) print("cluster probs p(c_m)") print(p_c) # update input means # TODO do this computation as a vector thing for m in range(M): input_means[m] = cw_expectation(x, m) print("input_means") print(input_means) # update input variances for m in range(M): input_vars[m] = cw_expectation((x - input_means)**2.0, m) + tiny print("input_vars") print(input_vars) # update functions for m in range(M): a = np.array([cw_expectation(y, m), cw_expectation(x * y, m)]) print("a") print(a) B = np.array([ [1.0, input_means[m]], [input_means[m], cw_expectation(x**2.0, m)] ]) print("B") print(B) beta[m] = np.linalg.inv(B).dot(a) print("beta") print(beta) # update output variances for m in range(M): predicted = np.array([cluster_func(x[n], m) for n in range(N)]) out_var[m] = cw_expectation((y - predicted)**2.0, m) + tiny print("output variances") print(out_var) for _ in range(10): em_step()
_psets/11.md 0 → 100644
 --- title: Problem Set 11 (Density Estimation) --- ## 1 {:.question} Generate a data set by sampling 100 points from \$\$y = \tanh(x)\$\$ between \$\$x = −5\$\$ to 5, adding random noise with a magnitude of 0.25 to \$\$y\$\$. Fit it with cluster-weighted modeling, using a local linear model for the clusters, \$\$y = a + bx\$\$. Plot the data and the resulting forecast \$\$\langle y | x \rangle\$\$, uncertainty \$\$\langle \sigma_y^2 | x \rangle\$\$, and cluster probabilities \$\$p(x|c_m) p(cm)\$\$ as the number of clusters is increased (starting with one), and animate their evolution during the EM iterations.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment