二次元的肥宅科学家们与宇宙物理生物计算机的异世界转生后创建的高级研究所的日常站点 - 代码-Code
SYSC 5108 W25 Assignment 3.pdfSYSC 5108 W25 ssignment 3William Wang101191168Q1From the Beta–Binomial conjugacy principle, a uniform $\theta \sim \text{Beta}(1,1)$ prior combined with the Binomial likelihood (with $x$ successes out of $n$ yields a posterior distribution$$ p(\theta \mid x) = \text{Beta}\bigl(x + 1,\; n - x + 1\bigr). $$Q2(a)Weighted sum = $0.20 \cdot 1 + (-0.10) \cdot 0.2 + 0.25 \cdot 0.5 + 0.05 \cdot 0.7 = 0.34$(b)Threshold ($\theta = 1$) output =$$ \begin{cases} 1 & \text{if } 0.34 \ge 1 \\ 0 & \text{otherwise} \end{cases} $$$$ $$So output = 0.(c)Logistic ($\varphi$) output =$$ \frac{1}{1 + e^{-0.34}} \approx 0.58. $$(d)ReLU output = $\max(0,\,0.34) = 0.34.$Q3(a) Logistic (sigmoid) activationNeuron 3$a_3 = \sigma(z_3)= \tfrac{1}{1 + e^{-(1.0 \times 0.7) + (0.6 \times 0.3) + 0.1 = 0.7 + 0.18 + 0.1}} = \tfrac{1}{1 + e^{-0.98}} \approx 0.73.$Neuron 4$a_4 = \sigma(z_4)= \tfrac{1}{1 + e^{-(0.9 \times 0.7) + (-0.8 \times 0.3) + 0.5 = 0.63 - 0.24 + 0.5 }} = \tfrac{1}{1 + e^{-0.89}} \approx 0.71.$Neuron 5$a_5 = \sigma(z_5) = \tfrac{1}{1 + e^{-(0.3 \times 0.73) + (0.2 \times 0.71) + 0.9 \approx 0.219 + 0.142 + 0.9}} = \tfrac{1}{1 + e^{-1.261}} \approx 0.78.$Neuron 5’s output $\approx 0.78$.(b) ReLU activationNeuron 3$a_3 = \max(0,\,0.98) = 0.98.$Neuron 4$a_4 = \max(0,\,0.89) = 0.89.$Neuron 5$a_5 = \max(0,\,1.372) = 1.372.$Hence, Neuron 5’s output $\approx 1.37$.(c) Filling in the network diagram (ReLU case)Weights and biases :$w_{3,1} = 1.0,\; w_{3,2} = 0.6,\; b_3 = 0.1$$w_{4,1} = 0.9,\; w_{4,2} = -0.8,\; b_4 = 0.5$$w_{5,3} = 0.3,\; w_{5,4} = 0.2,\; b_5 = 0.9$Neuron 3Weighted sum: $z_3 = 1.0 \times 0.7 + 0.6 \times 0.3 + 0.1 = 0.98$Activation: $a_3 = \mathrm{ReLU}(0.98) = 0.98$Neuron 4Weighted sum: $z_4 = 0.9 \times 0.7 + (-0.8) \times 0.3 + 0.5 = 0.89$Activation: $a_4 = \mathrm{ReLU}(0.89) = 0.89$Neuron 5Weighted sum: $z_5 = 0.3 \times a_3 + 0.2 \times a_4 + 0.9 = 1.372$Activation: $a_5 = \mathrm{ReLU}(1.372) = 1.372$Q4(a) Network Output$$ \begin{aligned} &\text{Hidden neurons:} \\ &a_2 = \max(0, 0.2 \times 0.5) = 0.1, \\ &a_3 = \max(0, 0.3 \times 0.5) = 0.15, \\ &a_4 = \max(0, 0.4 \times 0.5) = 0.2. \\ &\text{Output neuron:} \\ &y = a_5 = \max\bigl(0, (0.5 \times 0.1) + (0.6 \times 0.15) + (0.7 \times 0.2)\bigr) \\ &\quad = \max(0, 0.05 + 0.09 + 0.14) \\ &\quad = \max(0, 0.28) \\ &\quad = 0.28. \end{aligned} $$(b) Error Values ($\delta$) for Each Processing Neuron$$ E = \frac{1}{2} (t - y)^2 = \frac{1}{2} (0.9 - 0.28)^2 = \frac{1}{2} (0.62)^2 = 0.1922. $$(c) Gradients of the Error w.r.t. Each Weight$$ \delta_5 = (y - t)\cdot \mathbf{1}_{z_5>0} = (0.28 - 0.9)\times 1 = -0.62. $$For hidden neurons $i \in \{2,3,4\}$:$$ \delta_i = \delta_5 \times w_{5,i}, $$thus$$ \delta_2 = -0.62 \times 0.5 = -0.31, \quad \delta_3 = -0.62 \times 0.6 = -0.372, \quad \delta_4 = -0.62 \times 0.7 = -0.434. $$(d) Weight Updates with $\alpha=0.1$The gradient formulas are:$$ \frac{\partial E}{\partial w_{5,i}} = \delta_5\,a_i, \quad \frac{\partial E}{\partial w_{i,1}} = \delta_i \, x. $$Output layer:$$ \begin{aligned} \frac{\partial E}{\partial w_{5,2}} &= \delta_5 \, a_2 = (-0.62)\times0.1 = -0.062,\\ \frac{\partial E}{\partial w_{5,3}} &= (-0.62)\times0.15 = -0.093,\\ \frac{\partial E}{\partial w_{5,4}} &= (-0.62)\times0.2 = -0.124. \end{aligned} $$Hidden layer:$$ \begin{aligned} \frac{\partial E}{\partial w_{2,1}} &= \delta_2 \times x = (-0.31)\times0.5 = -0.155,\\ \frac{\partial E}{\partial w_{3,1}} &= (-0.372)\times0.5 = -0.186,\\ \frac{\partial E}{\partial w_{4,1}} &= (-0.434)\times0.5 = -0.217. \end{aligned} $$(e) Recalculation of the Error with Updated WeightsWith learning rate $\alpha=0.1$, each weight update is$$ w \leftarrow w - \alpha \,\frac{\partial E}{\partial w}. $$Output layer:$$ \begin{aligned} w_{5,2}^{\text{new}} &= 0.5 - 0.1\times(-0.062) = 0.5 + 0.0062 = 0.5062,\\ w_{5,3}^{\text{new}} &= 0.6 - 0.1\times(-0.093) = 0.6 + 0.0093 = 0.6093,\\ w_{5,4}^{\text{new}} &= 0.7 - 0.1\times(-0.124) = 0.7 + 0.0124 = 0.7124. \end{aligned} $$Hidden layer:$$ \begin{aligned} w_{2,1}^{\text{new}} &= 0.2 - 0.1\times(-0.155) = 0.2 + 0.0155 = 0.2155,\\ w_{3,1}^{\text{new}} &= 0.3 - 0.1\times(-0.186) = 0.3 + 0.0186 = 0.3186,\\ w_{4,1}^{\text{new}} &= 0.4 - 0.1\times(-0.217) = 0.4 + 0.0217 = 0.4217. \end{aligned} $$Checking the New Error(Forward Pass with Updated Weights):$$ \begin{aligned} a_2' &= \max(0, 0.2155\times0.5) = 0.10775,\\ a_3' &= \max(0, 0.3186\times0.5) = 0.1593,\\ a_4' &= \max(0, 0.4217\times0.5) = 0.21085,\\ y' &= a_5' = \max\bigl(0, 0.5062\times0.10775 + 0.6093\times0.1593 + 0.7124\times0.21085\bigr) \approx 0.302. \end{aligned} $$$$ E' = \frac{1}{2} (t - y')^2 = \frac{1}{2} (0.9 - 0.302)^2 = 0.1788 < 0.1922. $$Thus, the error has decreased, indicating successful training for this step.Q5a) What value will this network output?Define the input:$$ [1, 2, 3, 4, 5, 6] $$The first layer is a convolution with a window size of 3 and stride 1. The filter has weights:$w_1 = 1, \quad w_2 = 1, \quad w_3 = 1 $ and a bias: $w_0 = 0.75$Since this is a simple sum followed by a bias addition, we slide the window over the input, computing the convolution at each position.For the first position using ([1,2,3]):$$ z_1 = (1 + 2 + 3) + 0.75 = 6.75 $$$$ a_1 = \max(0, 6.75) = 6.75 $$For the second position using ([2,3,4]):$$ z_2 = (2 + 3 + 4) + 0.75 = 9.75 $$$$ a_2 = \max(0, 9.75) = 9.75 $$For the third position using ([3,4,5]):$$ z_3 = (3 + 4 + 5) + 0.75 = 12.75 $$$$ a_3 = \max(0, 12.75) = 12.75 $$For the fourth position using ([4,5,6]):$$ z_4 = (4 + 5 + 6) + 0.75 = 15.75 $$$$ a_4 = \max(0, 15.75) = 15.75 $$At this point, the feature map after the convolutional layer is:$$ [6.75, 9.75, 12.75, 15.75] $$Next, we apply max pooling with kernel size 2 and stride 2.The first pooling window takes:$$ \max(6.75, 9.75) = 9.75 $$The second pooling window takes:$$ \max(12.75, 15.75) = 15.75 $$So, after max pooling, the output is:$$ [9.75, 15.75] $$Now, these two values go into the fully connected layer. There is only one neuron, which has:Bias: $w_4 = 1 $Weights: $w_5 = 0.2, \quad w_6 = 0.5$The neuron computes:$$ z_7 = w_4 + w_5 \cdot 9.75 + w_6 \cdot 15.75 $$Substituting the values:$$ z_7 = 1 + (0.2 \times 9.75) + (0.5 \times 15.75) $$Calculating step-by-step:$$ 0.2 \times 9.75 = 1.95 $$$$ 0.5 \times 15.75 = 7.875 $$$$ z_7 = 1 + 1.95 + 7.875 = 10.825 $$Applying ReLU:$$ \max(0, 10.825) = 10.825 $$So, the final output of the network is:$$ \boxed{10.825} $$b) calculate the 𝛿 for each neuron in the network.The loss function is mean squared error:$$ E = \frac{1}{2} (t - y)^2 $$Taking the derivative with respect to ( N_7 ):$$ \delta_7 = \frac{\partial E}{\partial N_7} = (y - t) $$Substituting values:$$ \delta_7 = 10.825 - 1 = 9.825 $$$$ \delta_5 = \delta_7 \times w_5 = 9.825 \times 0.2 = 1.965 $$$$ \delta_6 = \delta_7 \times w_6 = 9.825 \times 0.5 = 4.9125 $$Neuron 5 pools $ (N_1, N_2) $, where the max was $ N_2 = 9.75 $.So only $ \delta_2 $ receives $ \delta_5 $, while $ \delta_1 = 0 $.$$ \delta_2 = \delta_5 = 1.965, \quad \delta_1 = 0 $$Neuron 6 pools $ (N_3, N_4) $, where the max was $ N_4 = 15.75 $.So only$ \delta_4 $ receives $ \delta_6 $, while $ \delta_3 = 0 $.$$ \delta_4 = \delta_6 = 4.9125, \quad \delta_3 = 0 $$Results:$$ \delta_1 = 0 $$$$ \delta_2 = 1.965 $$$$ \delta_3 = 0 $$$$ \delta_4 = 4.9125 $$$$ \delta_5 = 1.965 $$$$ \delta_6 = 4.9125 $$$$ \delta_7 = 9.825 $$c)$$ N_2 = \text{ReLU}(w_0 + w_1\cdot2 + w_2\cdot3 + w_3\cdot4). $$Taking partial derivatives:$$ \frac{\partial N_2}{\partial w_0} = 1, \quad \frac{\partial N_2}{\partial w_1} = 2, \quad \frac{\partial N_2}{\partial w_2} = 3, \quad \frac{\partial N_2}{\partial w_3} = 4. $$Since (N_2 > 0), the ReLU derivative is 1, so the gradient contribution from ( N_2 ) is:$$ \frac{\partial E}{\partial w_0} \Big|_{N_2} = \delta_2 \times 1 = 1.965, $$$$ \frac{\partial E}{\partial w_1} \Big|_{N_2} = \delta_2 \times 2 = 3.93, $$$$ \frac{\partial E}{\partial w_2} \Big|_{N_2} = \delta_2 \times 3 = 5.895, $$$$ \frac{\partial E}{\partial w_3} \Big|_{N_2} = \delta_2 \times 4 = 7.86. $$Since ( N_4 ) is computed as:$$ N_4 = \text{ReLU}(w_0 + w_1\cdot4 + w_2\cdot5 + w_3\cdot6). $$Taking partial derivatives:$$ \frac{\partial N_4}{\partial w_0} = 1, \quad \frac{\partial N_4}{\partial w_1} = 4, \quad \frac{\partial N_4}{\partial w_2} = 5, \quad \frac{\partial N_4}{\partial w_3} = 6. $$Since (N_4 > 0), the ReLU derivative is 1, so the gradient contribution from ( N_4 ) is:$$ \frac{\partial E}{\partial w_0} \Big|_{N_4} = \delta_4 \times 1 = 4.9125, $$$$ \frac{\partial E}{\partial w_1} \Big|_{N_4} = \delta_4 \times 4 = 19.65, $$$$ \frac{\partial E}{\partial w_2} \Big|_{N_4} = \delta_4 \times 5 = 24.5625, $$$$ \frac{\partial E}{\partial w_3} \Big|_{N_4} = \delta_4 \times 6 = 29.475. $$Since (w_j) appears in both (N_2) and (N_4), the total gradient for each weight is the sum of both contributions:$$ \frac{\partial E}{\partial w_0} = 1.965 + 4.9125 = 6.8775, $$$$ \frac{\partial E}{\partial w_1} = 3.93 + 19.65 = 23.58, $$$$ \frac{\partial E}{\partial w_2} = 5.895 + 24.5625 = 30.4575, $$$$ \frac{\partial E}{\partial w_3} = 7.86 + 29.475 = 37.335. $$Using the standard gradient descent update rule:$$ \Delta w_j = -\,\alpha \,\frac{\partial E}{\partial w_j}. $$So explicitly:$$ \Delta w_0 = -\,\alpha \times 6.8775, $$$$ \Delta w_1 = -\,\alpha \times 23.58, $$$$ \Delta w_2 = -\,\alpha \times 30.4575, $$$$ \Delta w_3 = -\,\alpha \times 37.335. $$Q6Feedforward (1 hidden layer, 30 nodes, KK outputs, input dimension DD):$\boxed{\text{Parameters} = (30D + 30) + (30K + K) = 30D + 30 + 31K.}$$$ (5 \times 5 \times 1) + 1 = 26 \text{ parameters}. $$$$ 26 \times 32 = 832. $$$$ (5 \times 5 \times 32) + 1 = 801. $$$$ 801 \times 64 = 51{,}264. $$$$ (3136 \times 1024) + 1024 = 3{,}212{,}288. $$$$ (1024 \times 10) + 10 = 10{,}250. $$$$ \begin{aligned} &\text{Conv1: } 832 \\ &\text{Conv2: } 51{,}264 \\ &\text{FC1: } 3{,}212{,}288 \\ &\text{FC2: } 10{,}250 \\ \hline &\textbf{Total: } 3{,}274{,}634. \end{aligned} $$Thus, the CNN has:$$ \boxed{3{,}274{,}634 \text{ parameters in total}.} $$Q7Task Complexity: Larger filters capture broader context but increase computation; smaller filters learn finer details.Input Image Size: Larger images may benefit from larger filters or bigger strides to reduce dimension quickly.Computational Constraints: More channels and bigger filters mean more parameters and higher memory/computational cost.Depth of the Network: Deeper networks often use smaller filters (e.g., 3×3) and fewer strides in early layers to preserve spatial information.
Q1Given n independent sample points $ x_i, i = 1, \dots, n $ from a Laplace distribution with PDF$$ f(x) = \frac{\lambda}{2} e^{-\lambda |x|} $$Provide MLE estimators for the parameter $ \lambda $.The log-likelihood for $n$ i.i.d. samples $x_1,\dots,x_n$ from the Laplace distribution$$ f(x)=\frac{\lambda}{2} e^{-\lambda|x|}\space \Rightarrow \space \ell(\lambda) =\sum_{i=1}^n \ln\left(\frac{\lambda}{2}\right) - \lambda\sum_{i=1}^n |x_i| = n \ln(\lambda) - n\ln(2) - \lambda \sum_{i=1}^n |x_i|. $$Differentiate and set to zero:$$ \frac{d\,\ell(\lambda)}{d\lambda} = \frac{n}{\lambda} - \sum_{i=1}^n |x_i| = 0 \quad\Longrightarrow\quad \hat{\lambda} = \frac{n}{\displaystyle \sum_{i=1}^n |x_i|}. $$The second derivative$\frac{d^2\,\ell(\lambda)}{d\lambda^2}=-\frac{n}{\lambda^2}<0$confirms a maximum.Thus, the MLE for $\lambda$ is$$ \boxed{\hat{\lambda} = \frac{n}{\sum_{i=1}^n |x_i|}}. $$Q2Maximum entropy is often used to select an a priori distribution when there is very little information about the RV. Consider the Bernoulli distribution $ p(x) = \begin{cases} p, & x = 1 \\ 1 - p, & x = 0 \\ 0, & \text{otherwise} \end{cases} $. Find its entropy $ H(X) $ and prove analytically that $ H(X) $ achieves its maximum when $ p = \frac{1}{2} $.For a Bernoulli random variable $X$ with $\Pr(X=1)=p$ and $\Pr(X=0)=1-p$, its entropy is$$ H(X)=-p\ln p-(1-p)\ln(1-p). $$Differentiate$$ \frac{d\,H}{dp} = -\ln p - 1 + \ln(1-p) + 1 = \ln\!\bigl(\tfrac{1-p}{p}\bigr) = 0 \quad\Longrightarrow\quad \tfrac{1-p}{p}=1 \quad\Longrightarrow\quad p=\tfrac12. $$Second derivative $ \tfrac{d^2\,H}{dp^2}=-\tfrac{1}{p} - \tfrac{1}{1-p} $ is negative for $0<p<1$, implying a maximum at $p=\tfrac12$.Q3A multivariate logistic regression model has been built to predict the propensity of shoppers to perform a repeat purchase of a free gift that they are given. The input features used by the model are the age of the customer, the socioeconomic band to which the customer belongs (a, b, or c), the average amount of money the customer spends on each visit to the shop, and the average number of visits the customer makes to the shop per week. This model is being used by the marketing department to determine who should be given the free gift. The weights in the trained model are shown in the following table. It is assumed that one hot coding is used for a, b, and c. Furthermore, if both b and c are 0, a will be used. Therefore, a is encoded as $b = c = 0$.FeatureWeightIntercept ($w[0]$)-3.82398AGE-0.02909Socioeconomic BAND B-0.00980Socioeconomic BAND C-0.19558SHOP VALUE0.02999SHOP FREQUENCY0.74572Use this model to make predictions for each of the following query instances.IDAGESocioeconomic BANDSHOP FREQUENCYSHOP VALUE156b1.60109.32221c4.9211.28348b1.21161.19437c0.72170.65532a1.08165.39import math # Given weights: w0 = -3.82398 w_age = -0.02909 w_band_b = -0.09089 w_band_c = -0.19558 w_value = 0.02999 w_frequency = 0.74572 # Helper function to compute sigmoid def sigmoid(z): return 1.0 / (1.0 + math.exp(-z)) # (Age, Band, ShopFrequency, ShopValue) queries = [ (56, 'b', 1.60, 109.32), (21, 'c', 4.92, 11.28), (48, 'a', 1.21, 161.19), (37, 'c', 0.72, 170.65), (32, 'a', 1.08, 165.39) ] # One-hot encoding for socioeconomic band: a => (0,0), b => (1,0), c => (0,1) def encode_band(band): if band == 'a': return (0,0) elif band == 'b': return (1,0) elif band == 'c': return (0,1) # Compute predictions for i, (age, band, freq, val) in enumerate(queries, start=1): b_enc, c_enc = encode_band(band) logit = (w0 + w_age * age + w_band_b * b_enc + w_band_c * c_enc + w_value * val + w_frequency * freq) # Probability prob = sigmoid(logit) print(f"Query {i}: logit = {logit:.4f}, predicted probability = {prob:.4f}")Query 1: logit = -1.0723, predicted probability = 0.2550Query 2: logit = -0.6232, predicted probability = 0.3490Query 3: logit = 0.5161, predicted probability = 0.6262Query 4: logit = 0.5588, predicted probability = 0.6362Query 5: logit = 1.0106, predicted probability = 0.7331Regression model:$$ \text{logit} = w_0 + w_{\text{age}}\cdot\text{Age} + w_{\text{bandB}}\cdot I_{\{\text{band=B}\}} + w_{\text{bandC}}\cdot I_{\{\text{band=C}\}} + w_{\text{value}}\cdot\text{ShopValue} + w_{\text{freq}}\cdot\text{ShopFrequency}. $$Then the predicted probability is the sigmoid of this linear combination:$$ = \frac{1}{1 + e^{-\text{logit}}}. $$For each row, substitute the features (with appropriate one‐hot encoding for the band) and compute for the five queries are:$p \approx 0.255$$p \approx 0.349$$p \approx 0.626$$p \approx 0.637$$p \approx 0.731$Q4In building multivariate logistic regression models, it is recommended that all continuous input features be normalized to the range [-1,1]. The following table shows a data quality report for the dataset used to train the model described in the last question.Data Quality ReportFeature%CountMiss.Card.Min.Qrt.MeanMedianQrt.Max.Std. Dev.AGE%5,200640182232.732326312.2SHOP FREQUENCY%5,20003160.21.02.21.31.34.35.4SHOP VALUE%5,20003,730511.8101.9100.14100.14174.6230.7Feature%CountMiss.Card.ModeMode Count%SOCIOECONOMIC BAND%5,20083a2,66451.2REPEAT PURCHASE%5,20002no2,79153.7On the basis of the information in this report, all continuous features were normalized using range normalization, and any missing values were replaced using mean for imputing continuous features and mode for imputing categorical features. After applying these data preparation operations, a multivariate logistic regression model was trained to give the weights shown in the following table.Model WeightsFeatureWeightIntercept ($w[0]$)0.6679AGE-0.5795SOCIOECONOMIC BAND B-0.1981SOCIOECONOMIC BAND C-0.2318SHOP VALUE3.4091SHOP FREQUENCY2.0499Use this model to make predictions for each of the query instances shown in the following table (question marks refer to missing values).Query InstancesIDAGESocioeconomic BANDSHOP FREQUENCYSHOP VALUE138a1.90165.39256b1.60109.32318c6.0010.094?b1.33204.62562?0.85110.50import numpy as np import pandas as pd # Given data min_age, max_age = 18, 63 min_shop_freq, max_shop_freq = 0.2, 4.3 min_shop_value, max_shop_value = 5, 174.6 mean_age = 32.7 mode_socioeconomic_band = "a" # Model weights intercept = 0.6679 weights = { "AGE": -0.5795, "SOCIOECONOMIC_BAND_B": -0.1981, "SOCIOECONOMIC_BAND_C": -0.2318, "SHOP_VALUE": 3.4091, "SHOP_FREQUENCY": 2.0499, } # Query data (before imputation and normalization) query_instances = pd.DataFrame([ [1, 38, 'a', 1.90, 165.39], [2, 56, 'b', 1.60, 109.32], [3, 18, 'c', 6.00, 10.09], [4, np.nan, 'b', 1.33, 204.62], [5, 62, np.nan, 0.85, 110.50], ], columns=["ID", "AGE", "SOCIOECONOMIC_BAND", "SHOP_FREQUENCY", "SHOP_VALUE"]) #missing values query_instances["AGE"].fillna(mean_age, inplace=True) query_instances["SOCIOECONOMIC_BAND"].fillna(mode_socioeconomic_band, inplace=True) # Normalize continuous features to [-1,1] range query_instances["AGE"] = 2 * (query_instances["AGE"] - min_age) / (max_age - min_age) - 1 query_instances["SHOP_FREQUENCY"] = 2 * (query_instances["SHOP_FREQUENCY"] - min_shop_freq) / (max_shop_freq - min_shop_freq) - 1 query_instances["SHOP_VALUE"] = 2 * (query_instances["SHOP_VALUE"] - min_shop_value) / (max_shop_value - min_shop_value) - 1 # Encode categorical variables query_instances["SOCIOECONOMIC_BAND_B"] = (query_instances["SOCIOECONOMIC_BAND"] == 'b').astype(int) query_instances["SOCIOECONOMIC_BAND_C"] = (query_instances["SOCIOECONOMIC_BAND"] == 'c').astype(int) # Compute logistic regression predictions def sigmoid(x): return 1 / (1 + np.exp(-x)) # Compute linear combination of inputs and weights linear_combination = ( intercept + weights["AGE"] * query_instances["AGE"] + weights["SOCIOECONOMIC_BAND_B"] * query_instances["SOCIOECONOMIC_BAND_B"] + weights["SOCIOECONOMIC_BAND_C"] * query_instances["SOCIOECONOMIC_BAND_C"] + weights["SHOP_FREQUENCY"] * query_instances["SHOP_FREQUENCY"] + weights["SHOP_VALUE"] * query_instances["SHOP_VALUE"] ) # Apply sigmoid function to get probabilities query_instances["PREDICTION"] = sigmoid(linear_combination) query_instances"ID", "PREDICTION"Query 1: $0.968$Query 2: $0.551$Query 3: $0.826$Query 4: $0.987$Query 5: $0.388$Q5Consider the training for logistic regression model. Assume $K = 2$, $y = 0$ or $1$. The loss function is defined below.$$ J(w) = -y \log(\sigma(z)) - (1 - y) \log(1 - \sigma(z)) $$where $z = w^T x, x_0 = 1$.(a) Show that$$ \nabla_w J(w) = -x(y - \sigma(z)) $$$$ z = w^T x, \quad \sigma(z) = \frac{1}{1+e^{-z}}, \quad J(w) = -\,y \,\log\!\bigl(\sigma(z)\bigr) \;-\; \bigl(1-y\bigr)\,\log\!\bigl(1-\sigma(z)\bigr). $$We want $\nabla_{w} J(w)$. A standard result (or by direct differentiation) shows:$$ \frac{\partial}{\partial z}\,J(w) = -\,\frac{y}{\sigma(z)}\,\sigma'(z) \;-\;\frac{1-y}{1-\sigma(z)}\,[-\sigma'(z)], $$where $\sigma'(z) = \sigma(z)\bigl(1 - \sigma(z)\bigr)$.Simplify:$$ \frac{\partial}{\partial z}\,J(w) = -\,y \,\frac{\sigma'(z)}{\sigma(z)} \;+\;(1-y)\,\frac{\sigma'(z)}{1-\sigma(z)}. $$$$ \frac{\partial}{\partial z}\,J(w) = \sigma(z)\bigl(1-\sigma(z)\bigr)\! \Bigl( -\,\frac{y}{\sigma(z)} \;+\; \frac{1-y}{\,1-\sigma(z}\Bigr). $$A little algebra shows that this combination inside the parentheses becomes $\sigma(z) - y$. Indeed:$$ -\,y \,\frac{1}{\sigma(z)} + (1-y)\,\frac{1}{1-\sigma(z)} = \frac{-\,y\,[1-\sigma(z)] + (1-y)\,\sigma(z)}{\sigma(z)\,\bigl(1-\sigma(z)\bigr)} = \frac{\sigma(z) - y}{\sigma(z)\,\bigl(1-\sigma(z)\bigr)}. \frac{\partial J}{\partial z} = \sigma(z) - y. $$Finally, by the chain rule with $z = w^T x$ and $\nabla_z z = x$:$$ \nabla_{w}\,J(w) \;=\; \bigl[\sigma(z) - y\bigr]\,x \;=\; -\,x\,(y \;-\;\sigma(z)). $$(b)Now change the encoding for $y$ to $y = -1$, or $1$. Show the result again. Explain the impact of this change.$$ \nabla_{\mathbf{w}} J(\mathbf{w}) = -\,\mathbf{x}\,[y - \sigma(z)]. $$Now, we switch the labels to $y \in \{-1, +1\}$. Below is the new loss function and the corresponding gradient derivation.Loss Function ($y \in \{-1, +1\}$)When $y$ takes values $\pm 1$, we can replace the original $\{0,1\}$ indicator with the following:$$ \mathbb{1}[\,y=+1\,] = \frac{1 + y}{2}, \quad \mathbb{1}[\,y=-1\,] = \frac{1 - y}{2}. $$$$ \boxed{ J(\mathbf{w}) = -\,\frac{1 + y}{2}\,\ln\bigl(\sigma(z)\bigr) - \frac{1 - y}{2}\,\ln\bigl(1 - \sigma(z)\bigr), } $$where $z = \mathbf{w}^\top \mathbf{x}$ and $\sigma(z) = \frac{1}{1+e^{-z}}$.Gradient Derivation Preparation$$ z = \mathbf{w}^\top \mathbf{x}, \quad \sigma(z) = \frac{1}{1 + e^{-z}}. $$To compute the gradient of the loss function $J(\mathbf{w})$ with respect to $\mathbf{w}$, we first compute the derivative of $z$, and then multiply by $\frac{\partial z}{\partial \mathbf{w}} = \mathbf{x}$:$$ \nabla_{\mathbf{w}} J(\mathbf{w}) = \frac{\partial J}{\partial z} \cdot \frac{\partial z}{\partial \mathbf{w}} = \left( \frac{\partial J}{\partial z} \right) \mathbf{x}. $$Derivation of $\frac{\partial J}{\partial z}$$$ \begin{aligned} J(\mathbf{w}) &= -\frac{1 + y}{2}\,\ln\bigl(\sigma(z)\bigr) - \frac{1 - y}{2}\,\ln\bigl(1-\sigma(z)\bigr), \\ \frac{\partial J}{\partial z} &= -\frac{1 + y}{2}\,\underbrace{\frac{\partial}{\partial z}\bigl[\ln(\sigma(z))\bigr]}_{\frac{1}{\sigma(z)}\,\sigma'(z)} - \frac{1 - y}{2}\,\underbrace{\frac{\partial}{\partial z}\bigl[\ln(1-\sigma(z))\bigr]}_{\frac{1}{1-\sigma(z)}\cdot[-\sigma'(z)]}. \end{aligned} $$$$ \sigma'(z) = \sigma(z)\,[1 - \sigma(z)], $$$$ \frac{\partial}{\partial z}\ln(\sigma(z)) = \frac{\sigma'(z)}{\sigma(z)} = \frac{\sigma(z)[1 - \sigma(z)]}{\sigma(z)} = 1 - \sigma(z), $$$$ \frac{\partial}{\partial z}\ln\bigl(1 - \sigma(z)\bigr) = \frac{-\,\sigma'(z)}{1 - \sigma(z)} = -\,\frac{\sigma(z)[1 - \sigma(z)]}{1 - \sigma(z)} = -\,\sigma(z). $$we can write the equation as:$$ \begin{aligned} \frac{\partial J}{\partial z} &= -\frac{1 + y}{2}\,\Bigl[\,1 - \sigma(z)\Bigr] + \frac{1 - y}{2}\,\bigl[-\,\sigma(z)\bigr] \\ &= -\frac{1 + y}{2}\;\Bigl[\,1 - \sigma(z)\Bigr] + \frac{1 - y}{2}\;\sigma(z). \end{aligned} $$$$ \begin{aligned} \frac{\partial J}{\partial z} &= -\frac{1 + y}{2} + \frac{1 + y}{2}\,\sigma(z) + \frac{1 - y}{2}\,\sigma(z) \\ &= -\frac{1 + y}{2} + \sigma(z)\,\Bigl[\,\frac{1 + y}{2} + \frac{1 - y}{2}\Bigr] \\ &= -\frac{1 + y}{2} + \sigma(z)\,\Bigl[\frac{(1 + y)+(1 - y)}{2}\Bigr] \\ &= -\frac{1 + y}{2} + \sigma(z)\,\Bigl[\frac{2}{2}\Bigr] \\ &= -\frac{1 + y}{2} + \sigma(z) \\ &= \boxed{\sigma(z) - \frac{1 + y}{2}}. \end{aligned} $$Substituting into the GradientUsing the chain rule:$$ \nabla_{\mathbf{w}} J(\mathbf{w}) = \left( \sigma(z) - \frac{1 + y}{2} \right) \mathbf{x}. $$Alternatively, we can write it in the "negative" form:$$ \nabla_{\mathbf{w}} J(\mathbf{w}) = -\,\mathbf{x}\,\Bigl[\frac{1 + y}{2} - \sigma(z)\Bigr]. $$Original labels $\{0,1\}$The gradient is:$$ \nabla_{\mathbf{w}} J(\mathbf{w}) = -\,\mathbf{x}\,[y - \sigma(z)]. $$New labels $\{-1,+1\}$The gradient is:$$ nabla_{\mathbf{w}} J(\mathbf{w}) = \left[ \sigma(z) - \frac{1 + y}{2} \right]\,\mathbf{x}. $$Here, $\frac{1 + y}{2}$ equals 1 when $y = +1$ and equals 0 when $y = -1$Q6Consider the training for logistic regression model using LSE. Assume $K = 2$, and the loss function is defined through the errors in probabilities as below.$$ J(\beta) = \sum_{i=1}^{M} \left[ I_1(y^{(i)}) \left( 1 - p_1(x^{(i)}) \right)^2 + I_2(y^{(i)}) \left( p_1(x^{(i)}) \right)^2 \right] $$where $z = \beta^\top x, x_0 = 1$. Find the gradient of the loss function and compare with the MLE discussed in class.$$ K = 2 \quad\Longrightarrow\quad \text{class 1 or class 2}. $$We define indicator functions$$ I_1(y^{(i)})= \begin{cases} 1 & \text{if }y^{(i)}=\text{class 1},\\ 0 & \text{otherwise}, \end{cases} \qquad I_2(y^{(i)})= \begin{cases} 1 & \text{if }y^{(i)}=\text{class 2},\\ 0 & \text{otherwise}. \end{cases} $$The model is logistic regression:We write$$ z^{(i)} = \beta^\top x^{(i)}, \quad p_1\bigl(x^{(i)}\bigr) = \sigma\bigl(z^{(i)}\bigr) = \frac{1}{1+e^{-z^{(i)}}}. $$Here $p_1(x^{(i)})$ is the predicted probability that $x^{(i)}$ belongs to class 1.In words:If the true label is class 1, the “error” term is $\bigl(1 - p_1(x^{(i)})\bigr)^2$.If the true label is class 2, the “error” term is $\bigl(p_1(x^{(i)})\bigr)^2$.Focus on the $i$-th data point. Its contribution to $J(\beta)$ is:$$ J_i(\beta) =\; I_1(y^{(i)})\;\bigl(1 - p_1^{(i)}\bigr)^2 \;+\; I_2(y^{(i)})\;\bigl(p_1^{(i)}\bigr)^2, $$where we set $p_1^{(i)} := p_1\bigl(x^{(i)}\bigr) = \sigma\bigl(\beta^\top x^{(i)}\bigr).$We will compute $\nabla_{\beta} J_i(\beta)$. By the chain rule:$$ \frac{\partial p_1^{(i)}}{\partial \beta} = \sigma'\bigl(\beta^\top x^{(i)}\bigr)\,x^{(i)} = p_1^{(i)}\bigl(1 - p_1^{(i)}\bigr)\,x^{(i)}. $$Hence,$$ \begin{aligned} \nabla_{\beta} \Bigl[\bigl(1 - p_1^{(i)}\bigr)^2\Bigr] &= 2\bigl(1 - p_1^{(i)}\bigr)\,\nabla_{\beta}\bigl(1 - p_1^{(i)}\bigr) \;=\; 2\bigl(1 - p_1^{(i)}\bigr)\,\bigl[-\,\nabla_{\beta} p_1^{(i)}\bigr] \\[4pt] &= -2\,\bigl(1 - p_1^{(i)}\bigr)\,p_1^{(i)}\bigl(1 - p_1^{(i)}\bigr)\,x^{(i)} \;=\; -2\,p_1^{(i)}\,\bigl(1 - p_1^{(i)}\bigr)^2\,x^{(i)}. \\ \nabla_{\beta} \Bigl[\bigl(p_1^{(i)}\bigr)^2\Bigr] &= 2\,p_1^{(i)}\,\nabla_{\beta} p_1^{(i)} = 2\,p_1^{(i)}\,p_1^{(i)}\bigl(1 - p_1^{(i)}\bigr)\,x^{(i)} = 2\,\bigl(p_1^{(i)}\bigr)^2\,\bigl(1 - p_1^{(i)}\bigr)\,x^{(i)}. \end{aligned} $$Therefore,$$ \nabla_{\beta} J_i(\beta) = I_1\bigl(y^{(i)}\bigr) \bigl[-2\,p_1^{(i)}\,\bigl(1 - p_1^{(i)}\bigr)^2\,x^{(i)}\bigr] \;+\; I_2\bigl(y^{(i)}\bigr) \bigl[\,2\,\bigl(p_1^{(i)}\bigr)^2\,\bigl(1 - p_1^{(i)}\bigr)\,x^{(i)}\bigr]. $$Factor out $2\,p_1^{(i)}\bigl(1 - p_1^{(i)}\bigr)\,x^{(i)}$:$$ \begin{aligned} \nabla_{\beta} J_i(\beta) &= 2\,p_1^{(i)}\bigl(1 - p_1^{(i)}\bigr)\,x^{(i)} \Bigl[ I_2\bigl(y^{(i)}\bigr)\,p_1^{(i)} \;-\; I_1\bigl(y^{(i)}\bigr)\,\bigl(1 - p_1^{(i)}\bigr) \Bigr] \\[6pt] &= 2\,p_1^{(i)}\bigl(1 - p_1^{(i)}\bigr)\,x^{(i)} \Bigl[ -\,I_1\bigl(y^{(i)}\bigr)\,\bigl(1 - p_1^{(i)}\bigr) \;+\; I_2\bigl(y^{(i)}\bigr)\,p_1^{(i)} \Bigr]. \end{aligned} $$Now notice:If $y^{(i)}=\text{class 1}$, then $I_1(y^{(i)})=1$, $I_2(y^{(i)})=0$. The bracket becomes:$$ -\,\bigl(1 - p_1^{(i)}\bigr) \;=\; p_1^{(i)} \;-\;1. $$If $y^{(i)}=\text{class 2}$, then $I_1(y^{(i)})=0$, $I_2(y^{(i)})=1$. The bracket becomes:$$ p_1^{(i)}. $$In fact, that bracket is exactly$$ p_1^{(i)} - I_1\bigl(y^{(i)}\bigr) $$in both cases. Therefore, a more compact final expression is:$$ \boxed{ \nabla_{\beta} J_i(\beta) \;=\; 2\,p_1^{(i)}\bigl(1 - p_1^{(i)}\bigr)\, \Bigl[ p_1^{(i)} \;-\; I_1\bigl(y^{(i)}\bigr) \Bigr] \,x^{(i)}. } $$Summing over all $i=1,\dots,M$,$$ \boxed{ \nabla_{\beta} J(\beta) \;=\; \sum_{i=1}^M \nabla_{\beta} J_i(\beta) \;=\; 2\,\sum_{i=1}^M p_1^{(i)}\bigl(1 - p_1^{(i)}\bigr)\,\Bigl[p_1^{(i)} - I_1\bigl(y^{(i)}\bigr)\Bigr]\,x^{(i)}. } $$where $p_1^{(i)} := \sigma\bigl(\beta^\top x^{(i)}\bigr)$.Comparison with MLE (Cross-Entropy) GradientFor logistic regression with the standard cross-entropy loss,$$ \text{(Cross-Entropy)} \quad J_{\text{CE}}(\beta) \;=\; -\,\sum_{i=1}^M \Bigl[ I_1\bigl(y^{(i)}\bigr)\,\ln\bigl(p_1^{(i)}\bigr) \;+\; I_2\bigl(y^{(i)}\bigr)\,\ln\bigl(1 - p_1^{(i)})\bigr) \Bigr], $$the gradient is famously$$ \boxed{ \nabla_{\beta} J_{\text{CE}}(\beta) = \sum_{i=1}^M \Bigl[ p_1^{(i)} - I_1\bigl(y^{(i)}\bigr) \Bigr]\, x^{(i)}. } $$ConclusionGradient for LSE Loss$$ \nabla_{\beta} J(\beta) = 2 \sum_{i=1}^M p_1^{(i)}\bigl(1 - p_1^{(i)}\bigr) \Bigl[ p_1^{(i)} \;-\; I_1(y^{(i)}) \Bigr] x^{(i)}. $$Gradient for MLE (Cross-Entropy) Loss$$ \nabla_{\beta} J_{\text{CE}}(\beta) = \sum_{i=1}^M \Bigl[ p_1^{(i)} \;-\; I_1(y^{(i)}) \Bigr] x^{(i)}. $$Q7The following images are handwritten instances of the digits 0 and 1. The images are small, 8 pixels by 8 pixels, and each pixel contains a gray level from the range [0, 7].Rather than use individual pixel values, which can lead to very high-dimensional feature vectors, a simpler way to represent images for use with regression models is to calculate a histogram for each image and use this as the feature vector instead. In this case, the histograms simply count the frequency of occurrence of each possible gray level in each image. The table that follows shows the histograms for a small dataset of 10 images split between examples of digits 0 and 1.IDGL-0GL-1GL-2GL-3GL-4GL-5GL-6GL-7DIGIT0313627564013731413213123134187370338230115141431532525110532632151516313523621207314341551118384222448193832344491A logistic regression model has been trained to classify digits as either 0 or 1. The weights in this model are as follows:InterceptGL-0GL-1GL-2GL-3GL-4GL-5GL-6GL-7w[0]w[1]w[2]w[3]w[4]w[5]w[6]w[7]w[8]0.3090.100-0.152-0.1630.191-0.631-0.716-0.478-0.171This model has been used to make predictions for the instances in the training set above. These predictions, and the related calculations required for calculating error and $w[j]$ error values are shown in the following table.IDσ(zi)yierrors.e.w[0]w[1]w[2]w[3]w[4]w[5]w[6]w[7]w[8]00.00010-0.0001?0.00000.00000.00000.00000.00000.00000.00000.00000.00001?10.85860.73730.10423.85620.31270.10420.41690.10420.31270.20841.354920.000000.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000?30.056510.94350.8902?1.91180.10060.15090.00000.05030.05030.25160.70444?0-0.00180.0000-0.00010.00000.00000.00000.00000.00000.00000.00000.000050.02561?0.94950.02430.77650.14560.07280.04850.02430.02430.12130.33976?0-0.00130.00000.00000.00000.00000.00000.00000.00000.00000.00000.000070.00450?0.00000.0000-0.0006?-0.0001-0.0001-0.0000-0.0001-0.0001-0.000280.020910.9791?0.02000.75980.08000.08000.04000.04000.08000.80000.160090.029210.97080.94250.02751.04470.08250.05500.08250.11000.11000.02750.2474In the table, the $w[j]$ errors are calculated as $\frac{\partial L(w)}{\partial w[j]} = (\gamma_i - \sigma(z_i))\sigma(z_i)(1-\sigma(z_i))x_i[j]$ which is the gradient for $W[l]$ estimated using least squared error (LSE) instead of MLE. $x_i[j]$ is the $j^{th}$ feature of the $i^{th}$ sample point.a)Some of the model predictions are missing in the preceding table (marked with a ?). Calculate these.ID = 1$ y_1 = 1 $$ \text{error}_1 = 0.8586 $Hence,$$ \sigma(z_1) = y_1 - \text{error}_1 = 1 - 0.8586 = 0.1414. $$ID = 4$ y_4 = 0 $$ \text{error}_4 = -0.0018 $Hence,$$ \sigma(z_4) = y_4 - \text{error}_4 = 0 - (-0.0018) = 0.0018. $$ID = 6$ y_6 = 0 $$ \text{error}_6 = -0.0013 $Hence,$$ \sigma(z_6) = y_6 - \text{error}_6 = 0 - (-0.0013) = 0.0013. $$Therefore, the three missing model predictions ($\sigma(z_i)$ values) are:ID 1: $\sigma(z_1) = 0.1414$ID 4: $\sigma(z_4) = 0.0018$ID 6: $\sigma(z_6) = 0.0013$.b)Some of the Error and Squared Error values are missing in the preceding table (marked with a ?). Calculate these values.ID = 0Given:$\sigma(z_0) = 0.0001$$y_0 = 0$$\text{error}_0 = -0.0001$ (already provided in the table)$\text{s.e.}_0 = ?$Thus:$$ \text{s.e.}_0 \;=\; (-0.0001)^2 \;=\; 0.00000001 \;\approx\; 0.0000 \quad \text{(rounded to four decimal places)} $$Hence the missing s.e. for ID=0 is 0.0000.ID = 5Given:$\sigma(z_5) = 0.0256$$y_5 = 1$$\text{error}_0 = ?$ (already provided in the table)$\text{s.e.}_0 = 0.9495$Thus:$$ \text{error}_5 \;=\; y_5 - \sigma(z_5) \;=\; 1 - 0.0256 \;=\; 0.9744. $$Hence the missing error for ID=5 is 0.9744.ID = 7Given:$\sigma(z_7) = 0.0045$$y_7 = 0$$\text{error}_7 = ?$$\text{s.e.}_7 = ?$$$ \text{error}_7 \;=\; y_7 - \sigma(z_7) \;=\; 0 - 0.0045 \;=\; -0.0045. $$Hence the missing values for ID=7 are error = -0.0045.ID = 8Given:$\sigma(z_8) = 0.0209$$y_8 = 1$$\text{error}_8 = 0.9791$ (already in the table)$\text{s.e.}_8 = ?$$$ \text{s.e.}_8 \;=\; (0.9791)^2 \;=\; 0.95864881 \;\approx\; 0.9586 \quad \text{(rounded to four decimal places).} $$Hence the missing s.e. for ID=8 is 0.9586.Final Answers the missing values:ID 0:$\text{s.e.} = 0.0000$ID 5:$\text{error} = 0.9744$ID 7:$\text{error} = -0.0045$ID 8:$\text{s.e.} = 0.9586$c)Some of the $w[j]$ error values are missing in the preceding table (marked with a ?). Calculate these.The $w[j]$ error for sample $i$:$$ \frac{\partial L(w)}{\partial w[j]} = (\gamma_i - \sigma(z_i)) \;\sigma(z_i)\bigl(1 - \sigma(z_i)\bigr)\; x_i[j]. $$$ \gamma_i $ is the true label.$ \sigma(z_i) $ is the model prediction.$ x_i[j] $ is the feature value.$ w[0] $ is the intercept ($ x_i[0] = 1 $).ID = 2, $w[8]$ error$ \sigma(z_2) = 0.0000 $, $ \gamma_2 = 0$$ \gamma_2 - \sigma(z_2) = 0 $$ \sigma(z_2)(1 - \sigma(z_2)) = 0 $$$ (\gamma_2 - \sigma(z_2))\,\sigma(z_2)(1-\sigma(z_2))\,x_2[8] = 0.0 \times 0.0 \times (\dots) = 0.0. $$$$ \boxed{w[8]\text{ error for ID=2} = 0.0000.} $$ID = 3, $w[0]$ error$ \sigma(z_3) = 0.0565 $, $ \gamma_3 = 1 $$ \gamma_3 - \sigma(z_3) = 0.9435 $$ \sigma(z_3)(1 - \sigma(z_3)) \approx 0.0533 $$$ \frac{\partial L}{\partial w[0]} = (0.9435) \times (0.0533) \times (1) \approx 0.0503. $$$$ \boxed{w[0]\text{ error for ID=3} = 0.0503.} $$ID = 7, $w[2]$ error$ \sigma(z_7) = 0.0045 $, $ \gamma_7 = 0 $$ \gamma_7 - \sigma(z_7) = -0.0045 $$ \sigma(z_7)(1-\sigma(z_7)) \approx 0.00448 $$$ (\gamma_7 - \sigma(z_7))\;\sigma(z_7)(1-\sigma(z_7)) = (-0.0045) \times 0.00448 \approx -0.00002016. $$$$ \frac{\partial L}{\partial w[0]} = -0.00002016 \times 1 = -0.00002016. $$Rounded: -0.0000.$$ \boxed{w[0]\text{ error for ID=7} = -0.0000.} $$Final Missing (w[j]) ErrorsID 2, (w[8]) = 0.0000ID 3, (w[0]) = 0.0503ID 7, (w[2]) = -0.0000d)Calculate a new set of weights for this model using a learning rate of 0.01 with LSE.New $w[j]$ errorIDσ(zi)yierrors.e.w[0]w[1]w[2]w[3]w[4]w[5]w[6]w[7]w[8]00.00010-0.00010.00000.00000.00000.00000.00000.00000.00000.00000.00000.000010.141410.85860.73730.10423.85620.31270.10420.41690.10420.31270.20841.354920.000000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.000030.056510.94350.89020.05031.91180.10060.15090.00000.05030.05030.25160.704440.00180-0.00180.0000-0.00010.00000.00000.00000.00000.00000.00000.00000.000050.025610.97440.94950.02430.77650.14560.07280.04850.02430.02430.12130.339760.00130-0.00130.00000.00000.00000.00000.00000.00000.00000.00000.00000.000070.00450-0.00450.00000.0000-0.0006-0.0000-0.0001-0.0001-0.0000-0.0001-0.0001-0.000280.020910.97910.95860.02000.75980.08000.08000.04000.04000.08000.80000.160090.029210.97080.94250.02751.04470.08250.05500.08250.11000.11000.02750.2474A single GD update with learning rate $ \alpha=0.01 $:$$ w[j] \;\leftarrow\; w[j] \;-\; \alpha\,G_j \;=\; w[j] \;-\; 0.01 \times G_j. $$$w[0]$:$$ 0.309 - 0.01 \times 0.2262 = 0.309 - 0.002262 \approx 0.3067. $$$w[1]$:$$ 0.100 - 0.01 \times 8.3484 = 0.100 - 0.083484 \approx 0.0165. $$$w[2]$:$$ -0.152 - 0.01 \times 0.7214 = -0.152 - 0.007214 \approx -0.1592. $$$w[3]$:$$ -0.163 - 0.01 \times 0.4628 = -0.163 - 0.004628 \approx -0.1676. $$$w[4]$:$$ 0.191 - 0.01 \times 0.5878 = 0.191 - 0.005878 \approx 0.1851. $$$w[5]$:$$ -0.631 - 0.01 \times 0.3288 = -0.631 - 0.003288 \approx -0.6343. $$$w[6]$:$$ -0.716 - 0.01 \times 0.5772 = -0.716 - 0.005772 \approx -0.7218. $$$w[7]$:$$ -0.478 - 0.01 \times 1.4087 = -0.478 - 0.014087 \approx -0.4921. $$$w[8]$:$$ -0.171 - 0.01 \times 2.8062 = -0.171 - 0.028062 \approx -0.1991. $$Final Updated Weights$w[0]$$w[1]$$w[2]$$w[3]$$w[4]$$w[5]$$w[6]$$w[7]$$w[8]$0.30670.0165-0.1592-0.16760.1851-0.6343-0.7218-0.4921-0.1991e)The following table shows handwritten examples of the digits 7 and 8 and their corresponding histogram values.IDGL-0GL-1GL-2GL-3GL-4GL-5GL-6GL-73515452483062054413iCalculate the output of the model (using the updated weights calculated in the previous part) for these two instances.$$ z \;=\; w[0]\cdot 1 \;+\; w[1]\cdot \text{GL-0} \;+\; \dots \;+\; w[8]\cdot \text{GL-7}, $$$$ \sigma(z) \;=\; \frac{1}{1 + e^{-z}}. $$For the digit “7”The histogram is:GL-0GL-1GL-2GL-3GL-4GL-5GL-6GL-7351545248Compute $z_7$:$$ \begin{aligned} z_7 &= w[0]\cdot 1 + w[1]\cdot 35 + w[2]\cdot 1 + w[3]\cdot 5 + w[4]\cdot 4 + w[5]\cdot 5 + w[6]\cdot 2 + w[7]\cdot 4 + w[8]\cdot 8 \\[6pt] &= 0.3067 + (0.0165)(35) + (-0.1592)(1) + (-0.1676)(5) + (0.1851)(4) + (-0.6343)(5) + (-0.7218)(2) + (-0.4921)(4) + (-0.1991)(8).\\ \end{aligned} \\ z_7 \approx -7.5489 \\ \sigma(z_7) = \frac{1}{1+e^{-(-7.5489)}} = \frac{1}{1+e^{7.5489}} \;\approx\; 0.00053. $$For the digit “8”The histogram is:GL-0GL-1GL-2GL-3GL-4GL-5GL-6GL-73062054413Compute $z_8$:$$ \begin{aligned} z_8 &= w[0]\cdot 1 + w[1]\cdot 30 + w[2]\cdot 6 + w[3]\cdot 2 + w[4]\cdot 0 + w[5]\cdot 5 + w[6]\cdot 4 + w[7]\cdot 4 + w[8]\cdot 13\\[6pt] &= 0.3067 + (0.0165)(30) + (-0.1592)(6) + (-0.1676)(2) + (0.1851)(0) + (-0.6343)(5) + (-0.7218)(4) + (-0.4921)(4) + (-0.1991)(13). \end{aligned} \\ z_8 \approx -11.1041. \\ \sigma(z_8) = \frac{1}{1 + e^{11.1041}} \;\approx\; 0.000015. $$iiComment on the appropriateness of these outputs.Using the updated logistic‐regression weights, both new images produce very negative (z)-values, hence very small sigmoid outputs:Predicted probability for “7” $\approx 0.00053$.Predicted probability for “8” $\approx 0.000015$.In other words, the model (after that single gradient‐descent step) strongly classifies both samples as closer to the “digit 0” class rather than “digit 1.”f)Now calculate new weights using MLE and one sample point ID = 0 only instead of step (d) and repeat step (e). Use the same learning rate.$$ \begin{aligned} \frac{\partial L}{\partial w[0]} &=\; 0.0001\times1 \,=\,0.0001,\\ \frac{\partial L}{\partial w[1]} &=\; 0.0001\times31 \,=\,0.0031,\\ \frac{\partial L}{\partial w[2]} &=\; 0.0001\times3 \,=\,0.0003,\\ \frac{\partial L}{\partial w[3]} &=\; 0.0001\times6 \,=\,0.0006,\\ \frac{\partial L}{\partial w[4]} &=\; 0.0001\times2 \,=\,0.0002,\\ \frac{\partial L}{\partial w[5]} &=\; 0.0001\times7 \,=\,0.0007,\\ \frac{\partial L}{\partial w[6]} &=\; 0.0001\times5 \,=\,0.0005,\\ \frac{\partial L}{\partial w[7]} &=\; 0.0001\times6 \,=\,0.0006,\\ \frac{\partial L}{\partial w[8]} &=\; 0.0001\times4 \,=\,0.0004. \end{aligned} $$$$ \begin{aligned} w[0] &\approx 0.309 - 0.01 \times 0.0001 = 0.309 - 0.000001 = 0.308999 \approx 0.3090,\\ w[1] &\approx 0.100 - 0.01 \times 0.0031 = 0.100 - 0.000031 = 0.099969 \approx 0.1000,\\ w[2] &\approx -0.152 - 0.01 \times 0.0003 = -0.152003 \approx -0.1520,\\ w[3] &\approx -0.163 - 0.01 \times 0.0006 = -0.163006 \approx -0.1630,\\ w[4] &\approx 0.191 - 0.01 \times 0.0002 = 0.190998 \approx 0.1910,\\ w[5] &\approx -0.631 - 0.01 \times 0.0007 = -0.631007 \approx -0.6310,\\ w[6] &\approx -0.716 - 0.01 \times 0.0005 = -0.716005 \approx -0.7160,\\ w[7] &\approx -0.478 - 0.01 \times 0.0006 = -0.478006 \approx -0.4780,\\ w[8] &\approx -0.171 - 0.01 \times 0.0004 = -0.171004 \approx -0.1710.\\ \end{aligned} $$Digit “7”Histogram: $\mathrm{GL{-}0}=35,\,\mathrm{GL{-}1}=1,\,\mathrm{GL{-}2}=5,\,\mathrm{GL{-}3}=4,\,\mathrm{GL{-}4}=5,\,\mathrm{GL{-}5}=2,\,\mathrm{GL{-}6}=4,\,\mathrm{GL{-}7}=8.$Using the (slightly) updated weights (w^{\text{new}}[j]) (to 6 decimals for accuracy), we get$$ \begin{aligned} z_{7} &\;=\;0.308999 \;+\; (0.099969)\times35 \;+\;(-0.152003)\times1 \;+\;(-0.163006)\times5 \;+\;(0.190998)\times4 \;+\;(-0.631007)\times5 \;+\;(-0.716005)\times2 \;+\;(-0.478006)\times4 \;+\;(-0.171004)\times8\\ &\;\approx\;-4.2621. \end{aligned} \\ \sigma(z_7) \;=\; \frac{1}{1+e^{-(-4.2621)}} \;=\; \frac{1}{1 + e^{4.2621}} \;\approx\; 0.014. $$Digit “8”Histogram: $\mathrm{GL{-}0}=30,\,\mathrm{GL{-}1}=6,\,\mathrm{GL{-}2}=2,\,\mathrm{GL{-}3}=0,\,\mathrm{GL{-}4}=5,\,\mathrm{GL{-}5}=4,\,\mathrm{GL{-}6}=4,\,\mathrm{GL{-}7}=13.$Similarly,$$ \begin{aligned} z_{8} &=\;0.308999 + (0.099969)\times30 + (-0.152003)\times6 + (-0.163006)\times2 + (0.190998)\times0 + (-0.631007)\times5 + (-0.716005)\times4 + (-0.478006)\times4 + (-0.171004)\times13\\ &\approx\;-8.0840. \end{aligned} \\ \sigma(z_8) \;=\; \frac{1}{1 + e^{8.0840}} \;\approx\;0.0003. $$ConclutionBecause the update from a single almost‐correct example $(y_0=0,\sigma(z_0)=0.0001)$ is vanishingly small,all new weights remain nearly the same as before (to four decimals).Consequently, the new model outputs for the “7” and “8” digits are still very close to zero—about 1.4% for the “7” and 0.03% for the “8.”Both are strongly classified as “digit 0” in this scenario
Assignment2.pdfQuestion 1: The Instrumental Variable ApproachStep 1) InitializationFirst, we define the number of time steps （$N$） along with the coefficients $a_1$ and $a_2$. To introduce noise, we set the standard deviation of white noise corresponding to variances of 4.0 and 9.0. Additionally, we initialize the arrays for the clean system output ($y$), the noisy system output ($y_m$), and the noise vector ($e$).% Simulation parameters N = 400; % Number of time steps a1 = 1.5; % Colored noise coefficient a2 = -0.7; % Colored noise coefficient sigma_e = 3; % Standard deviation of white noise (3 for variance = 9.0) rng(109); %set random seed % Initialize arrays y = zeros(N, 1); ym = zeros(N, 1); u = 2*(randi([0, 1], N, 1))-1; % Random binary input signal e = sigma_e * randn(N, 1); % White noise v = zeros(N, 1); % Initialize colored noiseStep 2) System SimulationFor the simulation, we first apply the clean filter:$$ y(k) = 1.5 y(k-1) - 0.7 y(k-2) + u(k-1) $$Then, we create the white noise vector:$$ e(k) = 2 \cdot \text{randn}(N,1) $$Finally, the actual system output is simulated:$$ y_m(k) = y(k) + e(k) $$We start the simulation from $k=3$ due to dependencies.Matlab Code% System simulation (start from k=3 due to dependencies) for k = 3:N y(k) = 1.5*y(k-1) - 0.7*y(k-2) + u(k-1); ym(k) = y(k) + e(k); % Measurement signal endStep 3) Initial Parameters for IV IterationsTo get the first instrumental solution, we compute the Least Squares (LS) estimate:$$ \theta_{LS} = (\varphi^T \varphi)^{-1} \varphi^T Y $$where $\varphi$ is the design matrix containing the regressors $[y_m(k-1), y_m(k-2), u(k-1)]$ and $Y$ is the output vector $y_m(3), y_m(4), ... y_m(N)$.Matlab Code% Design matrix for noise case Phi_noise_4 = [ym(2:end-1), ym(1:end-2), u(2:end-1)]; Y_noise_4 = ym(3:end); % Design matrix for clean case Phi_clean = [y(2:end-1), y(1:end-2), u(2:end-1)]; Y_clean = y(3:end);From the Matlab output, the LS estimations for clean, noisy (variance $4.0$), and noisy (variance $9.0$) cases are:CaseEstimated Parameters from LSClean Case (True Parameters)[1.500000, -0.700000, 1.000000]Noise Case (4.0 variance)[2.077900, -1.357300, 0.373200]Noise Case (9.0 variance)[2.949300, -2.518900, -0.316400]Step 4: Iterative IV EstimationAfter obtaining the initial estimates, the Instrumental Variable (IV) method is used to refine them iteratively. The instrumental variable is computed as:$$ X_i(k) = \theta_{LS}(1) X_i(k-1) - \theta_{LS}(2) X_i(k-2) + \theta_{LS}(3) u(k-1) $$At each iteration, the instrument matrix $Z$ is formed using lagged values of $X_i(k)$ and the input $u(k)$:$$ Z = [X_i(k-1), X_i(k-2), u(k-1)] $$The IV estimate is then updated as:$$ \theta_{\text{IV}} = (Z^T \Phi)^{-1} (Z^T Y) $$where $\Phi$ is the regressor matrix. The predicted output is:$$ y_{\text{pred}} = (\theta_{\text{IV}})^T \Phi_{\text{noise}} $$The cost function, measuring estimation error, is:$$ J = \sum (Y - y_{\text{pred}})^2 $$This iterative process improves parameter estimation by reducing the influence of noise.Matlab Code% Simulation parameters N = 400; a1 = 1.5; a2 = -0.7; sigma_e = 3; % Standard deviation of white noise (3 for variance = 9.0) rng(109); %set random seed % Initialize arrays y = zeros(N, 1); ym = zeros(N, 1); u = 2*(randi([0, 1], N, 1))-1; % Random binary input signal e = sigma_e * randn(N, 1); % White noise v = zeros(N, 1); % Initialize colored noise % System simulation (start from k=3 due to dependencies) for k = 3:N y(k) = 1.5*y(k-1) - 0.7*y(k-2) + u(k-1); ym(k) = y(k) + e(k); % Measurement signal end % Design matrix for noise case Phi_noise = [ym(2:end-1), ym(1:end-2), u(2:end-1)]; Y_noise = ym(3:end); % Design matrix for clean case Phi_clean = [y(2:end-1), y(1:end-2), u(2:end-1)]; Y_clean = y(3:end); % Initialize matrices to store the estimated parameters theta_IV_noise = zeros(3,5); theta_LS_noise = Phi_noise \ Y_noise; theta_IV_noise(:,1) = theta_LS_noise.'; cost = zeros(1,6); J_clean = zeros(1,6); %% Initialize plot part t = 3:N; % no noise nexttile; plot(t, y(t), 'LineWidth', 1.5); title('Noise-Free System Output y'); xlabel('Time Steps'); ylabel('y'); grid on; % With noise nexttile; plot(t, ym(t), 'LineWidth', 1.5); title('Noisy Measurement Signal ym'); xlabel('Time Steps'); ylabel('ym'); grid on; % create tiledlayout figure; tiledlayout(6,1, 'TileSpacing', 'Compact', 'Padding', 'Compact'); % Run IV iterations for iter = 1:10 if iter == 1 a1 = theta_LS_noise(1); a2 = theta_LS_noise(2); b1 = theta_LS_noise(3); else a1 = theta_IV_noise(1,iter-1); a2 = theta_IV_noise(2,iter-1); b1 = theta_IV_noise(3,iter-1); end Xi = zeros(N, 1); % Initialize instrumental variable for k = 3:N Xi(k) = a1 * Xi(k-1) + a2 * Xi(k-2) + b1 * u(k-1); end % Construct instrumental variable matrix Z = [Xi(2:end-1), Xi(1:end-2), u(2:end-1)]; % IV estimation theta_IV_noise(:,iter) = (Z' * Phi_noise) \ (Z' * Y_noise); % Predict output y_pred = Phi_clean * theta_IV_noise(:,iter); % Compute cost function cost(iter) = (1/(2*N)) * sum((Y_noise - y_pred).^2); J_clean(iter) = (1/(2*N)) * sum((Y_clean - y_pred).^2); % Display results fprintf('Iteration %d: Parameters: [%f, %f, %f], Cost: %f, J_clean: %f\n', ... iter, theta_IV_noise(:,iter), cost(iter),J_clean(iter)); nexttile; plot(t, y(t), 'k-', 'LineWidth', 0.1); hold on; plot(t, y_pred, 'k--', 'LineWidth', 0.1); title(sprintf('Iteration %d', iter)); xlabel('Time Steps'); ylabel('Output'); legend('No Noise y', sprintf('IV Predict y'), 'Location', 'best'); grid on; endPlot Variance = 4Plot Variance = 9ResultTrue Parameter [1.5, -0.7, 1]Estimated Parameters for Each Iteration (Noise Case: 4.0 Variance)Iteration$a_1$$a_2$$b_1$Cost$J$$J_{\text{clean}}$12.0779-1.35730.37322.28250.610021.4199-0.61440.54171.82680.110831.7092-0.90680.47141.87850.184841.5979-0.79410.49851.83770.135451.6385-0.83280.48891.84900.1501Estimated Parameters for Each Iteration (Noise Case: 9.0 Variance)Iteration$a_1$$a_2$$b_1$Cost$J$$J_{\text{clean}}$12.9493-2.5189-0.31647.78674.107720.9185-0.15240.49394.38650.447232.3535-1.6182-0.04575.09551.366541.3222-0.45500.35964.13090.264851.9760-1.14390.11044.39280.6198ConclusionInstrumental Variable (IV) estimation effectively reduces noise influence, but its performance depends on noise variance. Under moderate noise (4.0 variance), parameters converge smoothly toward the true values, with minor deviations. However, in high noise conditions (9.0 variance), convergence is slower, more erratic, and introduces bias, especially in $b_1$. The cost function decreases in both cases, but higher noise results in a consistently higher cost, indicating less accurate parameter estimates. Overall, IV estimation performs well under controlled noise but struggles with excessive noise.Question 2) Steepest Descent, LMS and RLS1. Problem SetupThe given system is represented in the Z-domain as:$$ \frac{d(z)}{x(z)} = \frac{0.2}{z - 0.8}. $$Rearranging:$$ zd(z) - 0.8d(z) = 0.2x(z). $$Converting to the time domain:$$ d(n + 1) = 0.2x(n) + 0.8d(n). $$Assume that the input signal $x(n)$ is white noise with variance $\text{var}[x(n)] = 1$.The parameter vector to estimate is:$$ \mathbf{H} = \begin{bmatrix} a \\ b \end{bmatrix}. $$2. Wiener--Hopf SolutionTo find the optimal parameters, we use the Wiener--Hopf equation:$$ \mathbf{H}_{\text{optimal}} = \mathbf{R}^{-1} \mathbf{p}. $$where:Cross-correlation vector:$$ \mathbf{p} = \mathbb{E}[ d(n) X(n) ], $$where $d(n)$ is the desired output, and $X(n)$ is the input vector.Autocorrelation matrix:$$ \mathbf{R} = \mathbb{E} [ X(n) X(n)^T ]. $$For this system:The input signals at time $n$ are $x(n)$ and $d(n)$, so the input vector is:$$ X(n) = \begin{bmatrix} d(n) \\ x(n) \end{bmatrix}. $$The cross-correlation vector is:$$ \mathbf{p} = \mathbb{E} \begin{bmatrix} d(n+1) d(n) \\ d(n+1) x(n) \end{bmatrix}. $$The autocorrelation matrix is:$$ \mathbf{R} = \begin{bmatrix} \mathbb{E}[ d^2(n) ] & \mathbb{E}[ d(n)x(n) ] \\[4pt] \mathbb{E}[ x(n)d(n) ] & \mathbb{E}[ x^2(n) ] \end{bmatrix}. $$With these expressions, we can simulate the system and precompute $\mathbf{p}$ and $\mathbf{R}$ for further calculations.Below is a straightforward way to carry out the Wiener--Hopf solution in closed form for the system$$ d(n+1) \;=\; 0.8\,d(n)\;+\;0.2\,x(n), $$with $x(n)$ white‐noise input of zero mean and variance $\sigma_x^2=1$. We seek the parameter vector$$ \mathbf{H} = \begin{bmatrix} a\\[3pt] b \end{bmatrix} $$that minimizes the mean‐square error between the "desired" $d(n+1)$ and the model output $\hat{d}(n+1) = a\,d(n)+b\,x(n)$.3. Autocorrelation Matrix $\mathbf{R}$We define the input vector$$ X(n) \;=\; \begin{bmatrix} d(n)\\[4pt] x(n) \end{bmatrix}. $$Then the autocorrelation matrix is$$ \mathbf{R} \;=\; \mathbb{E}\bigl[\,X(n)\,X(n)^\mathsf{T}\bigr] \;=\; \begin{bmatrix} \mathbb{E}[\,d(n)^2\,] & \mathbb{E}[\,d(n)\,x(n)\,]\\[3pt] \mathbb{E}[\,x(n)\,d(n)\,] & \mathbb{E}[\,x(n)^2\,] \end{bmatrix}. $$Compute $\mathbb{E}[x(n)^2]$.Since $x(n)$ is white noise with variance 1,$$ \mathbb{E}[\,x(n)^2\,] \;=\; 1. $$Compute $\mathbb{E}[\,d(n)\,x(n)\,]$.Because $d(n)$ depends on $x(n-1), x(n-2), \ldots$, and $x(n)$ is uncorrelated with past values of itself, it follows that $d(n)$ is uncorrelated with $x(n)$. Thus,$$ \mathbb{E}[\,d(n)\,x(n)\,] \;=\; 0. $$Compute $\mathbb{E}[\,d(n)^2\,]$.Denote $\sigma_d^2 = \mathbb{E}[\,d(n)^2\,]$. From the difference equation,$$ d(n+1) \;=\; 0.8 \, d(n) \;+\; 0.2 \, x(n). $$In steady state, using the fact that $d(n)$ is uncorrelated with $x(n)$ and$\sigma_x^2=1$,$$ \sigma_d^2 \;=\; \mathbb{E}[\,d(n)^2\,] \;=\; 0.8^2\,\sigma_d^2 \;+\; 0.2^2\,\underbrace{\mathbb{E}[\,x(n)^2\,]}_{1}. $$Hence$$ \sigma_d^2 \;=\; 0.64\,\sigma_d^2 \;+\; 0.04, \quad\Longrightarrow\quad 0.36\,\sigma_d^2 \;=\; 0.04 \quad\Longrightarrow\quad \sigma_d^2 \;=\; \frac{0.04}{0.36} \;=\; \tfrac{1}{9}. $$Therefore,$$ \mathbb{E}[\,d(n)^2\,] \;=\; \frac{1}{9}. $$Putting it all together,$$ \mathbf{R} \;=\; \begin{bmatrix} \frac{1}{9} & 0\\[4pt] 0 & 1 \end{bmatrix}. $$2. Cross‐Correlation Vector $\mathbf{p}$We need$$ \mathbf{p} \;=\; \mathbb{E}\bigl[\, d(n+1)\,X(n) \bigr] \;=\; \begin{bmatrix} \mathbb{E}[\,d(n+1)\,d(n)\,]\\[4pt] \mathbb{E}[\,d(n+1)\,x(n)\,] \end{bmatrix}. $$Use $d(n+1) = 0.8\,d(n) + 0.2\,x(n)$:First element $\mathbb{E}[\,d(n+1)\,d(n)\,]$:$$ \mathbb{E}[\,d(n+1)\,d(n)\,] \;=\; \mathbb{E}\bigl[\,(\,0.8\,d(n) + 0.2\,x(n)\,)\,d(n)\bigr] \;=\; 0.8\,\mathbb{E}[\,d(n)^2\,] \;+\; 0.2\,\mathbb{E}[\,x(n)\,d(n)\,]. $$Since $\mathbb{E}[\,x(n)\,d(n)\,] = 0$, we get$$ \mathbb{E}[\,d(n+1)\,d(n)\,] \;=\; 0.8 \times \frac{1}{9} \;=\; \frac{0.8}{9} \;=\; \frac{4}{45}. $$Second element $\mathbb{E}[\,d(n+1)\,x(n)\,]$:$$ \mathbb{E}[\,d(n+1)\,x(n)\,] \;=\; \mathbb{E}\bigl[\,(0.8\,d(n) + 0.2\,x(n))\,x(n)\bigr] \;=\; 0.8\,\underbrace{\mathbb{E}[\,d(n)\,x(n)\,]}_{0} \;+\; 0.2\,\mathbb{E}[\,x(n)^2\,] \;=\; 0.2 \times 1 \;=\; 0.2. $$Hence$$ \mathbf{p} \;=\; \begin{bmatrix} \frac{4}{45}\\[3pt] 0.2 \end{bmatrix}. $$4. Wiener--Hopf SolutionThe optimal parameter vector $\mathbf{H}_{\text{optimal}}$ is given by$$ \mathbf{H}_{\text{optimal}} \;=\; \mathbf{R}^{-1}\,\mathbf{p}. $$$$ \mathbf{R} \;=\; \begin{bmatrix} \frac{1}{9} & 0\\[3pt] 0 & 1 \end{bmatrix} \quad\Longrightarrow\quad \mathbf{R}^{-1} \;=\; \begin{bmatrix} 9 & 0\\[3pt] 0 & 1 \end{bmatrix}. $$Therefore,$$ \mathbf{H}_{\text{optimal}} \;=\; \begin{bmatrix} 9 & 0\\[3pt] 0 & 1 \end{bmatrix} \begin{bmatrix} \frac{4}{45}\\[3pt] 0.2 \end{bmatrix} \;=\; \begin{bmatrix} 9 \times \frac{4}{45}\\[3pt] 1 \times 0.2 \end{bmatrix} \;=\; \begin{bmatrix} \frac{36}{45}\\[3pt] 0.2 \end{bmatrix} \;=\; \begin{bmatrix} 0.8\\[3pt] 0.2 \end{bmatrix}. $$Hence the Wiener--Hopf solution for this system is$$ \boxed{ a \;=\; 0.8, \quad b \;=\; 0.2. } $$5.Steepest Descent AlgorithmMatlab Codeclear; clc; close all; %% Generate data N = 10000; % total samples x = randn(N,1); % white noise input, var=1 d = zeros(N,1); % desired signal initialization d(1) = 0; % initial condition % True system: d(n+1) = 0.8*d(n) + 0.2*x(n) for k = 1:N-1 d(k+1) = 0.8*d(k) + 0.2*x(k); end R = [1/9 0; 0 1 ]; p = [0.8/9; 0.2]; d2 = 1/9; H_opt = inv(R)*p; % -> [0.8 ; 0.2] %% Initialize Steepest Descent alpha = 0.01; % VERY TINY STEP size H_hat = [0; 0]; % Initial START [a; b] a_history = zeros(1,N); b_history = zeros(1,N); e_array = zeros(1,N); % store instantaneous error % Store initial values a_history(1) = H_hat(1); b_history(1) = H_hat(2); %% Steepest Descent Loop for n = 2:N e_array(n) = d(n) - H_hat(1)*d(n-1) - H_hat(2)*x(n-1); gradient = -2*p + 2*R*H_hat; H_hat = H_hat - alpha * gradient; a_history(n) = H_hat(1); b_history(n) = H_hat(2); end fprintf('Final estimated a = %.4f\n', a_history(end)); fprintf('Final estimated b = %.4f\n', b_history(end)); %% Contour Plot of J(a,b) + Parameter Path na = 60; % grid size for a-axis nb = 60; % grid size for b-axis a_range = linspace(-0.1,1.2,na); b_range = linspace(-0.1,0.4, nb); J_grid = zeros(na, nb); % Compute cost function grid for i = 1:na for j = 1:nb A = [a_range(i); b_range(j)]; J_grid(i,j) = d2 - 2*(A.'*p) + (A.'*R*A); end end figure; contour(a_range, b_range, J_grid.', 25); hold on; grid on; colorbar; xlabel('a'); ylabel('b'); title('Contour of J(a,b) + Steepest Descent Trajectory'); % Plot the parameter trajectory plot(a_history, b_history, 'r-','LineWidth',1.2,'MarkerSize',2); % Mark final estimate plot(a_history(end), b_history(end), 'ks','MarkerSize',6,'LineWidth',1.5); % Mark true Wiener solution plot(H_opt(1), H_opt(2), 'kx','MarkerSize',8,'LineWidth',2); legend('Cost contours','Path','Final estimate','True solution','Location','best'); hold off; %% plot a&b est figure figure; tiledlayout(2,1, 'TileSpacing', 'Compact', 'Padding', 'Compact'); % a_est plot nexttile; plot(a_history,'-','MarkerSize',2); hold on; yline(H_opt(1),'r--','LineWidth',1.5); xlabel('Iteration'); ylabel('a_{est}'); legend('a_{est}','true a=0.8','Location','best'); title('a_{est} vs. Iteration (Full)'); grid on; % b_est plot nexttile; plot(b_history,'-','MarkerSize',2); hold on; yline(H_opt(2),'r--','LineWidth',1.5); xlabel('Iteration'); ylabel('b_{est}'); legend('b_{est}','true b=0.2','Location','best'); title('b_{est} vs. Iteration (Full)'); grid on;6.LMS AlgorithmThe parameters are updated using the Least Mean Squares (LMS) rule:$$ e(n) = d(n) - [a_{\text{est}} d(n-1) + b_{\text{est}} x(n-1)] $$$$ H_{\text{hat}} = H_{\text{hat}} + \mu \cdot e(n) \cdot \begin{bmatrix} d(n-1) \\ x(n-1) \end{bmatrix} $$for n = 2:N % Input vector at time n (one-step-behind signals) Xn = [ d(n-1); x(n-1) ]; % Predicted output y_hat = H_hat.' * Xn; % = a_est*d(n-1) + b_est*x(n-1) % Error e_array(n) = d(n) - y_hat; % LMS parameter update H_hat = H_hat + mu * e_array(n) * Xn; % Store a_history(n) = H_hat(1); b_history(n) = H_hat(2); endPlotsWe will plot the following:Parameter Trajectory: $(a_{\text{est}}, b_{\text{est}})$Cost Surface Contour:The contour $J(a,b)$ derived from the Wiener formula.This allows visualization of how LMS moves within the cost surface.7.Standard Recursive Least Squares Algorithm(λ = 1.0)We use the classic RLS recursion:Input Vector:$$ X(n) = \begin{bmatrix} d(n-1) \\ x(n-1) \end{bmatrix} $$Predicted Output:$$ \hat{y}(n) = H_{\text{hat}}(n-1)^T X(n) $$Error Calculation:$$ e(n) = d(n) - \hat{y}(n) $$Gain Vector:$$ K(n) = \frac{P(n-1) X(n)}{\lambda + X(n)^T P(n-1) X(n)} $$Parameter Update:$$ H_{\text{hat}}(n) = H_{\text{hat}}(n-1) + K(n) e(n) $$Covariance Matrix Update:$$ P(n) = \frac{P(n-1) - K(n) X(n)^T P(n-1)}{\lambda} $$%% 4. RLS LOOP for n = 2:N % 1) Current regressor: X(n) = [ d(n-1); x(n-1) ] Xn = [ d(n-1); x(n-1) ]; % 2) Predicted output y_hat = H_hat.' * Xn; % 3) Error e_array(n) = d(n) - y_hat; % 4) Gain vector denom = lambda + Xn.' * P * Xn; % scalar K = (P * Xn) / denom; % 5) Parameter update H_hat = H_hat + K * e_array(n); % 6) Covariance update P = (P - K * (Xn.' * P)) / lambda; % Store results a_history(n) = H_hat(1); b_history(n) = H_hat(2); endPlotsWe visualize:The parameter trajectory $(a_{\text{est}}, b_{\text{est}})$.The cost surface contour $J(a, b)$ from the Wiener solution to compare the behavior of the RLS algorithm.8.Comment on the DifferenceSGD (Stochastic Gradient Descent)Performance: The core idea of SGD is to adjust the parameters by calculating the gradient of the error (that is, the direction in which the error changes the fastest), so that the parameters are updated in the direction of reducing the error. This method is very efficient in most cases, especially when the amount of data is large.Complexity: The calculation of SGD mainly relies on the multiplication of matrices and vectors. If the data scale is small (such as low dimensions of parameters and data), the amount of calculation is not large; but if the data scale is large (such as high-dimensional data or large-scale matrices), the calculation and storage will become very laborious. Therefore, although SGD itself is not complex, calculation and storage may become a bottleneck when processing large-scale data, especially when the signal is unstable (such as the data distribution changes over time).Implementation issues: SGD needs to store some data in advance (such as vector p and matrix R), which is feasible when the signal is stable. But if the signal is unstable (such as the data distribution changes over time) or the system characteristics change dynamically, the pre-stored data is no longer applicable.LMS (Least Mean Squares)Performance: LMS is a simpler gradient descent method, especially suitable for situations where data cannot be stored in advance (such as unstable signals or high-dimensional data). It does not directly calculate the entire error surface, but updates the parameters by estimating the current error, so it is called a "dirty gradient".Complexity: LMS is very simple to calculate because it only uses a fixed learning rate to update the parameters, avoiding complex matrix operations. This makes it very efficient and suitable for real-time processing.Implementation issues: The main problem with LMS is that it is very sensitive to the choice of learning rate. If the learning rate is set too small, the convergence speed will be slow; if it is set too large, it may cause parameter oscillation or even failure to converge.RLS (Recursive Least Squares)Performance: RLS is a very fast algorithm because it uses the inverse information of the autocorrelation matrix to adjust the parameters more accurately. It usually converges much faster than LMS and SGD.Complexity: The disadvantage of RLS is that it has a high computational complexity because it requires matrix operations, such as matrix inversion. These operations require high computing resources, especially when the data dimension is large.Implementation issues: Due to the high computational complexity of RLS, it is difficult to implement on hardware with limited resources. In addition, RLS has high requirements for numerical stability, and matrix inversion may cause numerical instability (such as "numerical expansion" phenomenon), which will affect the performance of the algorithm.Question 3a)Matlab Code% a) Recursive Least Squares (RLS) for system identification clear; clc; close all; rng(109); N = 400; sigma_e = 2; % Standard deviation of noise lambda_values = [1.0, 0.95]; % Forgetting factors theta_true = [1.5; -0.7; 1.0]; y = zeros(N,1); u = 2 * (rand(N,1) > 0.5) - 1; % Random binary input signal (±1) % Simulate the system for k = 3:N e = sigma_e * randn; y(k) = 1.5*y(k-1) - 0.7*y(k-2) + u(k-1) + e; end % Function to perform RLS estimation for lambda = lambda_values theta_est = zeros(3, N); P = 1000 * eye(3); trace_P = zeros(N,1); for k = 3:N phi = [y(k-1); y(k-2); u(k-1)]; % Regression vector K = P * phi / (lambda + phi' * P * phi); % Gain error = y(k) - phi' * theta_est(:, k-1); % Prediction error theta_est(:, k) = theta_est(:, k-1) + K * error; P = (P - (K * phi') * P) / lambda; % Update covariance trace_P(k) = trace(P); % Store trace end figure; subplot(2,1,1); plot(theta_est(1,:), 'r', 'LineWidth', 1.5); hold on; plot(theta_est(2,:), 'g', 'LineWidth', 1.5); plot(theta_est(3,:), 'b', 'LineWidth', 1.5); yline(1.5, 'r--', 'True \theta_1'); yline(-0.7, 'g--', 'True \theta_2'); yline(1.0, 'b--', 'True \theta_3'); xlabel('Time Step'); ylabel('\theta Estimates'); legend('\theta_1', '\theta_2', '\theta_3', 'Location', 'best'); title(['Parameter Estimates (\lambda = ', num2str(lambda), ')']); grid on; subplot(2,1,2); plot(trace_P, 'k', 'LineWidth', 1.5); xlabel('Time Step'); ylabel('Trace(P)'); title(['Trace of Covariance Matrix (\lambda = ', num2str(lambda), ')']); grid on; endIs there a difference?Yes!When comparing the results of $\lambda = 1.0$ and $\lambda = 0.95$, we find that they have significant differences in parameter stability and covariance matrix behavior. When $\lambda = 1.0$, the parameter estimates quickly stabilize and closely track the true values, with little fluctuation, and the trace of the covariance matrix also converges smoothly, indicating that the estimation process is stable and consistent. When $\lambda = 0.95$, the parameter estimates show more fluctuations because the algorithm pays more attention to recent data and is more sensitive to new observations. The trace of the covariance matrix also fluctuates more, reflecting the increased uncertainty in the estimation process.b)%% MATLAB Code for Part (b): RLS with Changing Input Signal clear all; close all; clc; rng(109); N = 400; sigma_e = 2; y = zeros(N,1); u = zeros(N,1); %% Generate the Input Signal u(k) % For k = 1 to 99, u(k) is random binary: +1 or -1 u(1:99) = 2*(rand(99,1) > 0.5) - 1; % For k = 100 to 300, set u(k) = +1.0 u(100:300) = 0; % For k = 301 to 400, u(k) returns to random binary: +1 or -1 u(301:400) = 2*(rand(100,1) > 0.5) - 1; %% Simulate the System y(1:2) = 0; for k = 3:N e = sigma_e * randn; y(k) = 1.5*y(k-1) - 0.7*y(k-2) + u(k-1) + e; end %% RLS Algorithm Setup lambda = 0.95; n_params = 3; theta_est = zeros(n_params, N); P = 1000 * eye(n_params); traceP = zeros(N,1); %% Run the RLS Algorithm for k = 3:N phi = [y(k-1); y(k-2); u(k-1)]; K = P * phi / (lambda + phi' * P * phi); error = y(k) - phi' * theta_est(:, k-1); theta_est(:, k) = theta_est(:, k-1) + K * error; P = (P - K * phi' * P) / lambda; traceP(k) = trace(P); end time = 1:N; % Time vector figure; subplot(2,1,1); plot(time, theta_est(1,:), 'r-', 'LineWidth', 2); hold on; plot(time, theta_est(2,:), 'b-', 'LineWidth', 2); plot(time, theta_est(3,:), 'g-', 'LineWidth', 2); xlabel('Time Step k'); ylabel('Parameter Estimates'); legend('\theta_1','\theta_2','\theta_3','Location','Best'); title('RLS Parameter Estimates with \lambda = 0.95'); grid on; subplot(2,1,2); plot(time, traceP, 'k-', 'LineWidth', 2); xlabel('Time Step k'); ylabel('Trace of Covariance Matrix P'); title('Trace of Covariance Matrix P over Time'); grid on; fprintf('Final parameter estimates:\n'); fprintf('theta1 = %.4f\n', theta_est(1,end)); fprintf('theta2 = %.4f\n', theta_est(2,end)); fprintf('theta3 = %.4f\n', theta_est(3,end));The RLS algorithm shows fast convergence, quickly stabilizing the parameter estimates ($\theta_1$, $\theta_2$, and $\theta_3$) . From chosen forgetting factor $\lambda = 0.95$. Between $k = 100$ and $k = 300$. The parameters remain stable as the input signal $u(k)$ is unchanged. However, at $k = 300$, when the input switches back to a binary signal, $\theta_3$ experiences significant fluctuations, indicating its higher sensitivity to input variations. The trace of the covariance matrix $P$ rises sharply between $k = 200$ and $k = 300$, reflecting increased uncertainty, but it starts decreasing after $k = 300$. It is showing that the algorithm is adapting to the new input conditions.c)%% MATLAB Code for Part (c): Time-Varying Parameters clear all; close all; clc; rng(109); N = 400; % Total number of data points sigma_e = 2; % Standard deviation of noise (sqrt(variance)=2) y = zeros(N,1); % Pre-allocate output vector u = zeros(N,1); % Pre-allocate input vector %% Generate the Input Signal u(k) % For k = 1 to 99, u(k) is random binary: +1 or -1 u(1:99) = 2*(rand(99,1) > 0.5) - 1; % For k = 100 to 300, set u(k) = +1.0 u(100:300) = 0; % For k = 301 to 400, u(k) returns to random binary: +1 or -1 u(301:400) = 2*(rand(100,1) > 0.5) - 1; %% Simulate the Time-Varying System y(1:2) = 0; for k = 3:N e = sigma_e * randn; if k < 250 % Before time step 250: original parameters y(k) = 1.5*y(k-1) - 0.7*y(k-2) + u(k-1) + e; else % At and after time step 250: changed parameters y(k) = 1.7*y(k-1) - 0.9*y(k-2) + 5*u(k-1) + e; end end %% RLS Algorithm Setup lambda = 0.95; n_params = 3; theta_est = zeros(n_params, N); P = 1000 * eye(n_params); traceP = zeros(N,1); %% Run the RLS Algorithm for k = 3:N phi = [y(k-1); y(k-2); u(k-1)]; K = P * phi / (lambda + phi' * P * phi); error = y(k) - phi' * theta_est(:, k-1); theta_est(:, k) = theta_est(:, k-1) + K * error; P = (P - K * phi' * P) / lambda; traceP(k) = trace(P); end %% Plot the Results time = 1:N; % Time vector figure; subplot(2,1,1); plot(time, theta_est(1,:), 'r-', 'LineWidth', 2); hold on; plot(time, theta_est(2,:), 'b-', 'LineWidth', 2); plot(time, theta_est(3,:), 'g-', 'LineWidth', 2); xlabel('Time Step k'); ylabel('Parameter Estimates'); legend('\theta_1','\theta_2','\theta_3','Location','Best'); title('RLS Parameter Estimates with a Parameter Change at k = 250'); grid on; subplot(2,1,2); plot(time, traceP, 'k-', 'LineWidth', 2); xlabel('Time Step k'); ylabel('Trace(P)'); title('Trace of the Covariance Matrix P over Time'); grid on; fprintf('Final parameter estimates:\n'); fprintf('theta1 = %.4f\n', theta_est(1,end)); fprintf('theta2 = %.4f\n', theta_est(2,end)); fprintf('theta3 = %.4f\n', theta_est(3,end));Observations on the Results:Parameter Estimates Before k = 250 :The estimated parameters ($\theta_1$, $\theta_2$, and $\theta_3$) are stable and remain relatively constant, showing that the RLS algorithm has successfully converged to an initial solution.$k$ = 250 :At $k$ = 250, the true system parameters were modified ($\theta_1$, $\theta_2$, and $\theta_3$). However, the estimated parameters take some time to reflect this change.Excitation at $k$ = 300At $k = 300$ the input signal $u(k)$ switches back to a random binary signal, introducing significant variations in the estimated parameters.$\theta_3$ (green) is particularly affected, experiencing a sharp increase, indicating that it is more sensitive to changes in the input signal.This suggests that $\theta_3$ is closely linked to the input excitation and that the sudden change in signal structure strongly influences it.Covariance Matrix Trace Behavior:The trace of the covariance matrix $P$ increases gradually from $k = 250$ and peaks around $k = 300$, which indicates a period of high uncertainty.Once the new parameter values are learned and the algorithm adapts, the trace begins to decrease, showing that confidence in the estimates is being restored.Key Information:The RLS algorithm successfully tracks parameter changes but takes time to adapt.The impact of the signal excitation at $k = 300$ is significant, especially for $\theta_3$.The covariance trace confirms increased uncertainty after parameter and input changes but stabilizes as the algorithm adapts.
import pandas as pd import os # change your path data_path = r'C:\Users\rtgtx\Desktop\SYSC5703\Destination_Prediction_Project' test_data = pd.read_csv(r'C:\Users\rtgtx\Desktop\SYSC5703\Destination_Prediction_Project\test.csv') submission_file_path = r'C:\Users\rtgtx\Desktop\SYSC5703\Destination_Prediction_Project\submission.csv' zone_feature_files_paths = { "Chukyo": os.path.join(data_path, "Zone_features", "Chukyo_zone_feature_area.csv"), "Higashisurugawan": os.path.join(data_path, "Zone_features", "Higashisurugawan_zone_feature_area.csv"), "Kinki": os.path.join(data_path, "Zone_features", "Kinki_zone_feature_area.csv"), "Kyushu": os.path.join(data_path, "Zone_features", "Kyushu_zone_feature_area.csv"), "Tokyo": os.path.join(data_path, "Zone_features", "Tokyo_zone_feature_area.csv") } train_files_paths = { "Chukyo": os.path.join(data_path, "train", "Chukyo.csv"), "Higashisurugawan": os.path.join(data_path, "train", "Higashisurugawan.csv"), "Kyushu": os.path.join(data_path, "train", "Kyushu.csv"), "Tokyo": os.path.join(data_path, "train", "Tokyo.csv") } # Load the data zone_feature_files = {name: pd.read_csv(path) for name, path in zone_feature_files_paths.items()} train_files = {name: pd.read_csv(path) for name, path in train_files_paths.items()} # Missing values def check_missing_values(dataframes): missing_values = {} for name, df in dataframes.items(): missing = df.isnull().sum() missing_values[name] = missing[missing > 0] return missing_values missing_zone_features = check_missing_values(zone_feature_files) missing_train_data = check_missing_values(train_files) missing_zone_features, missing_train_data import pandas as pd from sklearn.model_selection import train_test_split from sklearn.metrics import accuracy_score, confusion_matrix, classification_report from xgboost import XGBClassifier from sklearn.impute import SimpleImputer from sklearn.preprocessing import StandardScaler, LabelEncoder from sklearn.neighbors import KNeighborsClassifier def merge_zone_features(main_df, zone_df): return main_df.merge(zone_df, left_on='Origin', right_on='ZONE_ID', how='left') # Load and merge datasets merged_data = { "Chukyo": merge_zone_features(train_files['Chukyo'], zone_feature_files['Chukyo']), "Higashisurugawan": merge_zone_features(train_files['Higashisurugawan'], zone_feature_files['Higashisurugawan']), "Kyushu": merge_zone_features(train_files['Kyushu'], zone_feature_files['Kyushu']), "Tokyo": merge_zone_features(train_files['Tokyo'], zone_feature_files['Tokyo']) } combined_data = pd.concat(merged_data.values(), axis=0) # Top 100 most frequent destinations top_destinations = combined_data['Destination'].value_counts().head(100).index combined_data = combined_data[combined_data['Destination'].isin(top_destinations)] # Impute missing values using the mean strategy for numerical columns imputer = SimpleImputer(strategy='mean') combined_data'T000918002', 'T000918006', 'T000918021', 'T000918025', 'T000847001' = imputer.fit_transform( combined_data'T000918002', 'T000918006', 'T000918021', 'T000918025', 'T000847001' ) # Encode categorical variables label_encoders = {} for col in ['Gender', 'Occupation', 'Trip_type', 'Origin']: le = LabelEncoder() combined_data[col] = le.fit_transform(combined_data[col]) label_encoders[col] = le for col in ['Gender', 'Occupation', 'Trip_type']: # Get the labels seen during training known_labels = set(label_encoders[col].classes_) # Replace unseen labels with the most frequent label or a default value (e.g., -1 or 'unknown') test_data[col] = test_data[col].apply(lambda x: x if x in known_labels else np.nan) # Fill the NaN (unseen labels) with a default value (e.g., most frequent label seen during training) test_data[col].fillna(label_encoders[col].classes_[0], inplace=True) # Replace with most frequent label # Apply the label encoding test_data[col] = label_encoders[col].transform(test_data[col]) # Encode the target variable 'Destination' le_destination = LabelEncoder() y = le_destination.fit_transform(combined_data['Destination']) # Select features features = ['Gender', 'Age', 'Occupation', 'Trip_type', 'T000918002', 'T000918006', 'T000918021', 'T000918025', 'T000847001'] X = combined_data[features] # Scale features scaler = StandardScaler() X_scaled = scaler.fit_transform(X) print(combined_data['Destination'].nunique()) # Check number of unique values in 'Destination' print(combined_data['Destination'].value_counts().head(10)) # Inspect top 10 frequent classes # split X_train, X_test, y_train, y_test = train_test_split(X_scaled, y, test_size=0.2, random_state=64) # Train XGBoost with GPU xgb = XGBClassifier(use_label_encoder=False, eval_metric='mlogloss', tree_method='gpu_hist', random_state=64) xgb.fit(X_train, y_train) # Predict and evaluate y_pred = xgb.predict(X_test) accuracy = accuracy_score(y_test, y_pred) conf_matrix = confusion_matrix(y_test, y_pred) class_report = classification_report(y_test, y_pred) print("Accuracy:", accuracy) print("Confusion Matrix:\n", conf_matrix) print("Classification Report:\n", class_report) #skip this this is knn model X_train_knn, X_test_knn, y_train_knn, y_test_knn = train_test_split(X_scaled, y, test_size=0.2, random_state=48) knn = KNeighborsClassifier(n_neighbors=5) knn.fit(X_train, y_train) y_pred_knn = knn.predict(X_test_knn) accuracy = accuracy_score(y_test_knn, y_pred_knn) conf_matrix = confusion_matrix(y_test_knn, y_pred_knn) class_report = classification_report(y_test_knn, y_pred_knn) print("Training Accuracy:", accuracy) print("Confusion Matrix:\n", conf_matrix) print("Classification Report:\n", class_report) # SKip this from sklearn.model_selection import cross_val_score #XGBoost 5-fold cross-validation xgb_cv = XGBClassifier(use_label_encoder=False, eval_metric='mlogloss', tree_method='gpu_hist', random_state=48) scores = cross_val_score(xgb_cv, X_train, y_train, cv=5, scoring='accuracy') print("XGBoost Cross-Validation Accuracy Scores:", scores) print("Mean Cross-Validation Accuracy:", scores.mean()) import pandas as pd y_test_pred = xgb.predict(X_test_scaled) test_data['Destination'] = le_destination.inverse_transform(y_test_pred) submission = test_data'Trip_id', 'Destination' submission.to_csv(submission_file_path, index=False) print(f"Submission file saved to {submission_file_path}")