<jupyter_start><jupyter_text>**Read the `Instructions` notebook** before you start working on this problem set! It contains instructions on how to create the submission package and a detailed description of the provided classes (*BayesNet* and *Variable*).<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from bayesian_network import BayesNet<jupyter_output><empty_output><jupyter_text>Rejection SamplingImplement the rejection sampling algorithm and use it to approximate $P(A \mid D, E)$.In Rejection Sampling, one samples from the full joint distribution and throws away ('rejects') all samples that do not correspond to the evidence. ImplementationImplement- `sample_forward`, and- `rejection_sampling`.`sample_forward` returns a set of samples from the full joint distribution $P(\mathcal{X})$.`rejection_sampling` uses *sample_forward*, throws away all random events with mismatching evidence $\mathbf{E}$, and estimates the distribution of the query variable $X$. In other words, it computes an approximation of $P(X \mid \mathbf{E}=e)$. For simplicity, we will only consider probabilistic queries with one query variable $X$.We will need a function that samples from a discrete 1-D probability distribution. For this purpose, use the `sample_categorical()` function from the Instructions notebook, which takes a 1-D NumPy array representing a probability distribution and draws a sample from it:<jupyter_code>from utils import sample_categorical
help(sample_categorical)<jupyter_output>Help on function sample_categorical in module utils:

sample_categorical(dist: numpy.ndarray) -> numpy.int64
    Draws a single sample from a categorical distribution.
    :param dist: NumPy array listing the probability of each outcome.
    :returns: Index of the sampled element.<jupyter_text>Forward Sampling Implement the sample_forward function, which samples from the full joint distribution of a Bayesian network $\mathcal{B}$. (2 points)The `sample_forward` function must return one object:- samples from $\mathcal{B}$ of type `np.ndarray`, with shape `(sample_size, len(bayes_net))`**Hint**: Iterating over the *BayesNet* object returns the *Variable* objects in topological ordering. The probability distribution of variable $X$ given its parents $\mathit{pa}(X)$, $P(X \mid \mathit{pa}(X))$, can be obtained by passing the (possibly incomplete) random event to the variable, i.e., `variable(samples[i])`.<jupyter_code>def sample_forward(bayes_net: BayesNet, sample_size: int) -> np.ndarray:
    '''
    Samples from the full joint distribution.
    :param bayes_net: A Bayesian network of type BayesNet.
    :param sample_size: The number of samples to draw from the Bayesian network.
    :returns: A NumPy array of type np.int64 with shape (sample_size, len(bayes_net))
              containing samples from the Bayesian network
    '''
    
    # array holding the samples
    samples = np.empty((sample_size, len(bayes_net)), np.int64)
    
    # YOUR CODE HERE
    for i in range(sample_size):
        for variable in bayes_net:
            distribution = variable(samples[i])
            samples[i][variable.id] = sample_categorical(distribution)
    return samples
# sanity checks
bayes_net = BayesNet([(np.array([0.5, 0.5]), [0])])

samples = sample_forward(bayes_net, 3)

assert type(samples) == np.ndarray, f'\nWrong output type!\nExpected: np.ndarray\nGiven:\t  {type(samples)}'
assert samples.shape == (3, 1), f'\nWrong output shape!\nExpected: (3, 1)\nGiven:\t  {samples.shape}'
assert samples.dtype == np.int64, f'\nWrong numpy array data type!\nExpected: np.int64\nGiven:\t  {samples.dtype}'<jupyter_output><empty_output><jupyter_text>Rejection Sampling     Implement the rejection_sampling function, which estimates the probability distribution over the query variable $X$ given evidence  $\mathbf{E}=\mathbf{e}$, i.e. $P(X \mid \mathbf{E}=\mathbf{e})$. Use the sample_forward function implemented in the previous step to draw random events from $P(\mathcal{X})$. (3 points)The `rejection_sampling` function must return one object:-  The probability distribution over variable $X$ of type `np.ndarray` with shape `(bayes_net[query_variable].num_values,)`. If the sample from the FJD does not contain any random events with matching evidence, return a NumPy array with the same shape as $P(X \mid \mathbf{E}=\mathbf{e})$ containing only values `np.nan`, i.e. `np.full(bayes_net[query_variable].num_values, np.nan)`.**Hint**: Use `np.bincount` to count the number of occurrences of each value.<jupyter_code>def rejection_sampling(bayes_net: BayesNet, query_variable: int, evidence: dict={}, sample_size: int=100) -> np.ndarray:
    '''
    Estimates the distribution of the query variale given the value of the evidence variables.
    :param bayes_net: A Bayesian network of type BayesNet.
    :param query_variable: Id of the query variable (int).
    :param evidence: A dictionary of evidence variables (keys: int) and their correponding values (values: int).
    :param sample_size: The number of samples to use for the estimation.
    :returns: A NumPy array of type np.float64 representing the conditional distribution of the query variable given evidence, or a NumPy array of the same shape but with all values set to np.nan if no sample with matching evidence is found. 
    '''
    
    # forward sample...
    samples = sample_forward(bayes_net, sample_size)
    
    # YOUR CODE HERE
    for key, value in evidence.items():
        samples = samples[samples[:,key] == value]

    variable = bayes_net[query_variable]
    if len(samples) == 0:
        distribution = np.full(variable.num_values, np.nan)
    else:
        distribution = np.bincount(samples[:,query_variable], minlength=variable.num_values) / len(samples)
    
    return distribution

# sanity checks
bayes_net = BayesNet([(np.array([0.25, 0.25, 0.5]), [0])])

np.random.seed(0)
distribution = rejection_sampling(bayes_net, 0, {0:0}, 3)
assert distribution.shape == (3,), f'\nWrong output shape!\nExpected: (3,)\nGiven:\t  {distribution.shape}'
assert np.all(np.isnan(distribution))

distribution = rejection_sampling(bayes_net, 0, {}, 1000)
assert type(distribution) == np.ndarray, f'\nWrong output type!\nExpected: np.ndarray\nGiven:\t  {type(distribution)}'
assert distribution.shape == (3,), f'\nWrong output shape!\nExpected: (3,)\nGiven:\t  {distribution.shape}'
assert distribution.dtype == np.float64, f'\nWrong numpy array data type!\nExpected: np.float64\nGiven:\t  {distribution.dtype}'<jupyter_output><empty_output><jupyter_text>ExampleNow that we can approximate arbitrary probabilistic queries let us compare an approximation's accuracy to an exact solution.Consider the following Bayesian Network (all variables are binary):The conditional probability tables are given as:P(A)$a_0$$a_1$-0.20.8P(C)$c_0$$c_1$-0.90.1P(B | A)$a_0$$a_1$$b_0$0.90.2$b_1$0.10.8P(D | B, C)$b_0$$b_1$$c_0$$c_1$$c_0$$c_1$$d_0$0.10.20.990.8$d_1$0.90.80.010.2P(E | C)$c_0$$c_1$$e_0$0.70.4$e_1$0.30.6 First, let us store these tables as global NumPy arrays and create the BayesNet object:<jupyter_code>_A_, _B_, _C_, _D_, _E_ = 0, 1, 2, 3, 4

A = np.array([0.2, 0.8])
B_A = np.array([[0.9, 0.2], [0.1, 0.8]])
C = np.array([0.9, 0.1])
D_BC = np.array([[[0.1, 0.2], [0.99, 0.8]], [[0.9, 0.8], [0.01, 0.2]]])
E_C = np.array([[0.7, 0.4], [0.3, 0.6]])
              
bayes_net = BayesNet([
    (A, [_A_]),
    (B_A, [_B_, _A_]),
    (C, [_C_]),
    (D_BC, [_D_, _B_, _C_]),
    (E_C, [_E_, _C_])
])<jupyter_output><empty_output><jupyter_text>Exact Computation Compute $P(A \mid D, E)$ exactly (up to floating point precision). (1 point)Feel free to do this in a very inefficient manner, e.g., by first computing the joint probability distribution and then normalizing by the evidence.<jupyter_code>A_DE = None

# YOUR CODE HERE
fjdt = 1

for variable in bayes_net:
    fjdt = fjdt * variable.pdt

# Sum over the axis, that are not in our evidence
A_DE = fjdt.sum(axis=(_B_,_C_))
# Normalize by the sum of our A
A_DE = A_DE / A_DE.sum(axis=_A_)
# print("A_DE=\n", A_DE)
# print("A_DE.shape=", A_DE.shape)
# print("A_DE.sum(axis=0)=", A_DE.sum(axis=0))
# sanity checks
assert A_DE is not None
assert A_DE.shape == (2,2,2)
assert np.all(np.isclose(A_DE.sum(axis=0), 1))<jupyter_output><empty_output><jupyter_text>ComparisonRun the following code cell to plot the average [Kullback-Leibler divergence](https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence) between the exact distribution $P(A \mid D, E)$ and the approximations computed with rejection sampling. Different lines represent different value assignments to the evidence variables $D$ and $E$.**Hint**: The computation of the approximations might take a while...<jupyter_code>from utils import approx_error

np.random.seed(0)
errs = {}
sample_counts = np.array([10, 20, 40, 80, 160, 320, 640, 1280])

# compute approximation error
for d, e in zip([0, 0, 1, 1], [0, 1, 0, 1]):
    errs[(d, e)] = approx_error(
        bayes_net,
        rejection_sampling,  
        A_DE[:, d, e],
        _A_,
        {_D_:d, _E_:e}, 
        sample_counts,
        n_runs=100
    )

# plot
plt.figure(figsize=(10, 8))
plt.title('Kullback-Leiber Divergence')
for d, e in zip([0, 0, 1, 1], [0, 1, 0, 1]):
    plt.plot(sample_counts, errs[(d, e)], label=f'e:{e}, d:{d}', lw=2)
    
plt.legend()
plt.xscale('log')
plt.yscale('log')
plt.xlim(sample_counts.min(), sample_counts.max())
plt.xlabel('Number of samples')
plt.show()<jupyter_output><empty_output><jupyter_text>Answer the following question in one sentence! (1 point)Have a look at the average approximation errors when sampling $10^3$ samples with forward sampling. Why is the approximation error of $P(A \mid e=1, d=1)$ higher then the approximation error of $P(A \mid e=0, d=0)$? Because $P(e=1, d=1)$ is less likely than $P(e=0, d=0)$ and therefore it does not get sampled that often so the error is higher$ Store the answer to the following question into the provided result variable! (1 point)Assume $d=1$ and $e=0$. On average, how many of the random events sampled from $P(\mathcal{X})$ will be accepted by the rejection sampling algorithm (i.e., not rejected)? Give the fraction of samples accepted, e.g. 0.5, if every 2nd sample contains the correct evidence.<jupyter_code>result = None

# YOUR CODE HERE
p_d_e = fjdt.sum(axis=(_A_,_B_,_C_))
print(p_d_e)
p_d1_e0 = p_d_e[1,0]
print(p_d1_e0)
result = p_d1_e0
assert result is not None<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12261268143621085, 'fraction_numerical': 0.023968678380443087, 'mean_word_length': 4.726079825041006, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 35, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, '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': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '499261', 'license': 'Apache-2.0', 'n_tokens_mistral': 3568, 'n_tokens_neox': 3369, 'n_words': 1141, 'path': 'Problem 1.ipynb', 'repo_name': 'Dorkat0/Pommerman_students', 'size': '10472'}