<jupyter_start><jupyter_text>Density Estimation TechniquesIn this lab, you will implement the density estimation techniques we studied in the lecture.  Lab Instructions:1. Read the explanation above each requirement very well2. Read the requirement very well before jumping into the code.3. Some requirements have essay questions in them, make sure you do NOT miss them.4. PLEASE Read the hints! They are clear and made to help you complete the requirement as fast as you should About the LabThis lab consists of one requirement: [Requirement 1](Requirement-1) where you implement one of the density estimation methods.- Refer to the slides whenever you need an equation or a formula. - In this lab, you are free to choose your bump function. However, higher accuracies will get a higher grade. **You will need to justify your choice of that bump function, write that in your own words, do not copy from your colleagues**- Feel free to use code from Lab 3.  ------------------ Density Estimation RecapIn the previous lab, we were introduced to the problem of classification of multivariate data, and we applied Bayes classifier to classify these points. The Bayes classification rule depends mainly on the class-conditional probability (also known as the class density or class distribution). An accurate method to calculate this distribution is crucial for improving the performance of the classifier. We built our classifier based on the assumption that this data is multivariate Gaussian random distribution, and we estimated the mean and covariance of this distribution. Finally, a test vector x is classified as (i.e. assigned) the class that maximizes the posterior probability:$P(C_i|x_{test}) = \frac{P(x_{test}|C_i) * P(C_i)}{P(x_{test})}$Classification of $x_{test}$ = $argmax_{i} P\big(C_i|x_{test}\big)$ = $argmax_{i} P(x|C_i) * P(C_i)$* $P(C_i|x_{test})$: Posterior probability* $P(x|C_i)$: Class-conditional probability (or distribution)* $P(C_i)$: Class apriori probability                **Note that** $P(C_i|x_{test}) \neq P(x_{test}|C_i) * P(C_i)$. Instead,  $P(C_i|x_{test}) = \frac{P(x_{test}|C_i) * P(C_i)}{P(x_{test})}$. However, when we compare multiple classes, the denominator $P(x_{test})$ is independent of the class $i$ and can be regarded as normalizing factor.Where $P(x|C_i)$ is given by the multivariate Gaussian random distribution:$P(x|C_i) = \mathcal{N}(x; \mu_i, \Sigma_i) = \frac{1}{(2\pi)^{\frac{N}{2}}|\Sigma_i|^{\frac{1}{2}}} exp\big(\frac{-1}{2}(x-\mu_i)^T\Sigma^{-1}_{i}(x-\mu_i)\big)$  In this labWe will drop off the assumption that $P(x|C_i)$ is drawn out of a Gaussian distribution with mean $\mu$ and covariance $\Sigma$. The distribution is not necessarily Gaussian, yet we will still apply Bayes classifier to classify between data points by assigning each point to the class that maximizes its posterior probability as usual. What will change in Bayes classifier under the new assumption?We will answer that for you! $P(x|C_i)$ is no longer equal to the Gaussian distribution as shownin the previous equation i.e. we can’t use the Gaussian formula to substitute for $P(x|C_i)$. How will we figure out the unknown distribution (also known as likelihood function)?We will answer that for you too! There are three ways for doing so:    1. Histogram Analysis.    2. Naïve Estimator.    3. Parzen Window Density Estimator.    All these techniques aim to estimate a density function $P(x|C_i)$ for the training points, and hence you can use the density function to apply Bayes classifier and classify your test points.*Refer to the slides to get a brief about these techniques.*  Requirement 1You are given files **data.csv**, **test_data.csv**, and  **test_data_true.csv**. **The first file** contains a list of one-dimensional data points and their corresponding classes. The data points were drawn of the distribution shown in Figure 1 (NOTE: It will be your task to estimate this distribution). ![Screen%20Shot%202020-11-21%20at%2001.01.52-2.png](attachment:Screen%20Shot%202020-11-21%20at%2001.01.52-2.png)As shown in Figure 1, the points belonging to class 1 were drawn out of the distribution in blue,whereas the points belonging to class 2 were drawn out of the distribution in red.After applying Parzen Window Density Estimator and plotting the estimated distribution forboth classes, you should observe a similar distribution for each class to the correspondingdistribution shown in Figure 1.The format of data files is shown in the following table:||Class|Value||-|-|-||Point1| 1| 5.34||Point2| 2| 2.73||Point3| 2| -1.832|**The second file** contains test points without the associated classes. The true values for testpoints are found in *test_data_true.csv*.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import csv<jupyter_output><empty_output><jupyter_text>**The first step** is to read and import the training data, test data and the true values for the test data,found in the given files respectively.<jupyter_code>def read_data():
    data = None
    # TODO 1: Read the file 'data.csv' into the variable data.
    #        Convert the matrix data into numpy array and make sure that the values in data is float and not strings.
    # Hint: You already implemented it in Lab 3.
    # Hint: use data = np.asarray(data)
    # Note: data contains both the training data together with labelled classes.
    data = np.genfromtxt('data.csv', delimiter=',', dtype=float)
    return data

def read_test_data():
    # TODO 2: Read the file 'test_data.csv' and 'test_data_true.csv' into the variables test_data and test_data_true.
    #         Convert the matrix data into numpy arrays and make sure that the values in data is float and not strings.
    test_data = np.genfromtxt('test_data.csv', delimiter=',', dtype=float)

    # Hint: You already implemented it in Lab 3.
    # Note: test_data contains the unlabelled test class.
    # Note: test_data_true contains the actual classes of the test instances, which you will compare
    #       against your predicted classes.
    test_data_true = np.genfromtxt('test_data_true.csv', delimiter=',', dtype=float)
    return test_data, test_data_true

data = read_data()
test_data, test_data_true = read_test_data()<jupyter_output><empty_output><jupyter_text>**The second step** is to choose a suitable bump function (Parzen window), and find the optimal value of `h`for the bump function for each class.First, fill in the function `calculate_h_optimal()`. Pass the convenient parameters.<jupyter_code># Note: Pass the suitable parameters
def calculate_h_optimal(sigmai, N, M):
    h = 0
    # TODO 3: Calculate the optimal value for the bandwidth (h) of Parzen window (or the bump function you are going
    #  to compute). Pass the suitable parameters. 
    # Hint: Refer to the slides to find out the equation of computing h optimal and know which parameters you should pass.
    h = sigmai * (4 / ((2 + N) * M))**(1 / (N + 4))
    return h<jupyter_output><empty_output><jupyter_text>Second, implement the bump function (i.e. Parzen window). You should computethe bump function for any x, given a point p and an optimal value for h.Y = F(x), where x represents the domain of the function or equivalently somevalues to compute the bump function at, y represents the output of the bumpfunction and F represents the bump function.**Notes:** You should vectorise this (i.e. your function should work if x were a single data point or a vector of M data points).<jupyter_code>def bump_function(point, x, h):
    f = 0
    # TODO 4: Implement the bump function phi(x) at point p with bandwidth h. Note that x can be a scalar or a vector.
    # NOTE: Your solution should be vectorized.
    f = np.exp(-((x - point)**2) / (2 * h**2)) / (np.sqrt(2 * np.pi) * h)
    # f = 1/h * 1/((2*np.pi)**(1/2)) * np.exp(-np.dot((x-point), (x-point)) / (2 * h**2))
    return f<jupyter_output><empty_output><jupyter_text>**Question 1** What are the criteria of choosing bump function?<jupyter_code>'''
    Your Answer:
'''<jupyter_output><empty_output><jupyter_text>**Question 2** Justify why you choose that as your bump function. Do you expect it to enhance accuracy? Why? Write everything you would like to say<jupyter_code>'''
    Your Answer:
'''<jupyter_output><empty_output><jupyter_text>**At this point,** you need to iterate over the training points and apply the bumpfunction for every training point.\The density estimate will be the accumulation (i.e. the sum) of the output of these bump functions.We first need to generate the x vector which represents the possible values for the onedimensional feature we have. You can find the minimum and maximum possible values for x among the training data, and generate a range of points which will represent your one-dimensional axis.                            x = np.linspace(min_x, max_x, number_of_axis_points)                            Here, `min_x` represents the minimum value for the feature x, `max_x` represents themaximum value for the feature x, and `number_of_axis_points` represents theresolution of the scale (or how many points should be included between these twolimits).<jupyter_code># TODO 5: Construct x-axis that will be used in density estimation.
min_x = min(min(data[:, 1]), min(test_data)) - 5
max_x = max(max(data[:, 1]), max(test_data)) + 5
number_of_axis_points = 10000
x = np.linspace(min_x, max_x, number_of_axis_points)
print((x))<jupyter_output>[-12.02452718 -12.01972425 -12.01492132 ...  35.99039413  35.99519707
  36.        ]<jupyter_text>Now we need to compute the desnsities for each class. Remember that the density estimate represents `P(x|Ci)`We start by computing h_optimal for each class we have:<jupyter_code># TODO 6: Calculate the optimal h for each class. 
# (Warning: Do not use for loops. Use vectorization, broadcasting and masking).
h1 = calculate_h_optimal(np.std(data[data[:, 0] == 1.0, 1]), data.shape[1] - 1, data.shape[0])
h2 = calculate_h_optimal(np.std(data[data[:, 0] == 2.0, 1]), data.shape[1] - 1, data.shape[0])<jupyter_output><empty_output><jupyter_text>We then use these h_optimals to calculate the density functions:<jupyter_code># TODO 7: Calculate the probability density function p(x|C1) and p(x|C2) for each value on the constructed x-axis
pxgivenC1 = np.sum(np.array([bump_function(point, x, h1) for point in data[data[:, 0] == 1.0, 1]]), axis = 0) / len(data[data[:, 0] == 1.0]) # class 1
pxgivenC2 = np.sum(np.array([bump_function(point, x, h2) for point in data[data[:, 0] == 2.0, 1]]), axis = 0) / len(data[data[:, 0] == 2.0]) # class 2<jupyter_output>6<jupyter_text>Now we need to visualize how well are we estimating the densities.\Let's plot the density estimates for the two classes that we computed in **ONE PLOT** with **DIFFERENTCOLORS.** This is an ordinary plot.**Check your plot result**:\Is it similar to Figure 1? If yes, proceed to question 5. If no, unfortunately you havedone something wrong. Think: Do you need another bump function? Do you need more points? **FAQ:** What do you mean by two classes in one plot with different colors?\Answer: The two classes are overlaid together in one plot (such as Figure 1) withdifferent colors (red and blue). We can’t be any more specific :")<jupyter_code># plot the estimated densities
# TODO 8: Plot the two densities P(x|C1) and P(x|C2) on the same plot. use plt.plot(..)
plt.plot(x, pxgivenC1, label = 'P(x|C1)')
plt.plot(x, pxgivenC2, label = 'P(x|C2)')<jupyter_output><empty_output><jupyter_text>It’s time to classify the test points. You know what to do already! Given the points inthe test file, it’s required to classify each point whether it belongs to class 1 or 2 byapplying Bayes classifier. Report the accuracy of your classifier by comparing withthe true values found in the third file.Hint: It’s not as easy as it looks like, there’s a trick here. Be creative in your solution.Solutions may differ in this part, and this may lead to varying accuracies betweenstudents.\*May the odds be in your favor and you get the highest accuracy!*<jupyter_code># TODO 9: Calculate the aprior probabilities p(C1) and p(C2).
pc1 = len(data[data[:, 0] == 1.0]) / len(data)
pc2 = len(data[data[:, 0] == 2.0]) / len(data)
# TODO 10: Calculate the posterior probabilities p(C1|x) and p(C2|x) for all values in x.
pc1GivenX = pxgivenC1 * pc1 
pc2GivenX = pxgivenC2 * pc2

def find_nearest(array, value):
    '''
    This function finds which element in array is equal to value. If not found, it returns the nearest (or closest)
    element in the array to value.
    :param array: An array of values to search in.
    :param value: A value to select from the array.
    :return: The index of the found element in the array. If not found, it returns the index of the closest element.
    '''
    ## DO NOT Change Code ##
    array = np.asarray(array)
    index = (np.abs(array - value)).argmin()
    return index
# TODO 11: Compute the accuracy of your classifier.
accuracy = 0
for idx, point in enumerate(test_data):
    classification = 1 if pc1GivenX[find_nearest(x, point)] > pc2GivenX[find_nearest(x, point)] else 2
    if classification == test_data_true[idx]:
        accuracy += 1
    print("for point:", point, "the classification: ", classification)
accuracy /= len(test_data)

print("Accuracy = ", accuracy, "%")<jupyter_output>Accuracy =  0.91 %
{'fraction_non_alphanumeric': 0.09686609686609686, 'fraction_numerical': 0.01949317738791423, 'mean_word_length': 4.45562372188139, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 43, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '182505', 'license': 'MIT', 'n_tokens_mistral': 4072, 'n_tokens_neox': 3851, 'n_words': 1746, 'path': 'Lab 4 - Density Estimation/Density_Estimation.ipynb', 'repo_name': 'MuhammeedAlaa/Pattarn-Labs', 'size': '13338'}