How to determine the size of biases? I'm new to the world of machine learning. My question is how can I determine the size of the biases in a neural network (with backpropagation algorithm)? Currently, I have a 2-layer neural network (1 hidden and 1 output layer). Here's the code:

import numpy as np
from matplotlib import pyplot as plt 

sigmoid = lambda x : 1 / (1 + np.exp(-x))
dsigmoid = lambda y: y * (1 - sigmoid(y))

# This function performs the given function (func) to the whole numpy array
def mapFunc(array, func) :
    newArray = array.copy()
    for element in np.nditer(newArray, op_flags=['readwrite']) :
        element[...] = func(element)
    return newArray

class NeuralNetwork :

def __init__(self, input_nodes, hidden_nodes, output_nodes) :
    self.input_nodes = input_nodes
    self.hidden_nodes = hidden_nodes
    self.output_nodes = output_nodes

    self.W_ih = np.random.rand(hidden_nodes, input_nodes)
    self.W_ho = np.random.rand(output_nodes, hidden_nodes)

    self.B_ih = np.random.rand(hidden_nodes, 1)
    self.B_ho = np.random.rand(output_nodes, 1)

    self.learningRate = 0.1

def predict(self, inputs) :
    # Calculate hidden's output
    H_output = np.dot(self.W_ih, inputs)
    H_output += self.B_ih
    H_output = mapFunc(H_output, sigmoid) # Activation

    # Calculate output's output
    O_output = np.dot(self.W_ho, H_output)
    O_output += self.B_ho
    O_output = mapFunc(O_output, sigmoid) # Activation

    return O_output

def train(self, inputs, target) :
    # Calculate hidden's output
    H_output = np.dot(self.W_ih, inputs)
    H_output += self.B_ih
    H_output = mapFunc(H_output, sigmoid) # Activation

    # Calculate output's output
    O_output = np.dot(self.W_ho, H_output)
    O_output += self.B_ho
    O_output = mapFunc(O_output, sigmoid) # Activation

    # Calculate output error :
    O_error = O_output - target

    # Calculate output delta
    O_gradient = mapFunc(O_output, dsigmoid)
    O_gradient = np.dot(O_gradient, np.transpose(O_error)) * self.learningRate

    W_ho_delta = np.dot(O_gradient, np.transpose(H_output))

    self.W_ho -= W_ho_delta
    self.B_ho -= O_gradient

    # Calculate hidden error :
    W_ho_t = np.transpose(self.W_ho)
    H_error = np.dot(W_ho_t, O_error)

    # Calculate hidden delta :
    H_gradient = mapFunc(H_output, dsigmoid)
    H_gradient = np.dot(H_gradient, np.transpose(H_error)) * self.learningRate

    W_ih_delta = np.dot(H_gradient, inputs)

    self.W_ih -= W_ih_delta
    self.B_ih += H_gradient

    return O_output


n = NeuralNetwork(2, 2, 1)

inputs = np.matrix([[1], [0], [1], [1], [0], [1], [0], [0]])

input_list = []
input_list.append([[1], [0]])
input_list.append([[0], [1]])
input_list.append([[1], [1]])
input_list.append([[0], [0]])

target = np.matrix([[0], [0], [1], [1]])

outputs = []
for i in range(50000) :
    ind = np.random.randint(len(input_list))
    inp = input_list[ind]
    out = n.train(inp, target[ind]).tolist()
    outputs.append(out[0][0])

print outputs
plt.plot(outputs)

plt.show()

newInput = [[1], [1]]
print (n.predict(newInput))


In the train function, the line self.B_ih += H_gradient throws me an error about their sizes not being equal. I even tried to make the biases only a single number but that didn't help as it gets changed by H_gradient to a matrix. So, is there something wrong in the bias itself or I did some other step(s) wrong?
The biases appear in the activation function and can be seen as node feature. Suppose you have act = sigmoid, X your input layer and W your first hidden layer weight matrix.
Then the output matrix of this hidden layer will be act(X.dot(W)+b). 

So your bias should be a vector containing a random value for each node of the current hidden layer. 

It will shift randomly the weight and so the activation function. See here (usually it prevent from overfitting).

I didn't read the whole code so I don't know if the problem comes from here, but why don't you learn how to use tensorflow instead of re-write everything ? There are a lot of tutorials in the net, like MNIST.