neural network for binary classification of xor gate
There are a few mistakes in the code, so I am going to present a revised version here with comments.

Setup

import numpy as np
import pandas as pd
x=np.array([[0,0],[0,1],[1,0],[1,1]])
y=np.array([[0],[1],[1],[0]])
np.random.seed(0)

# Optional, but a good idea to have +ve and -ve weights
theta1=np.random.rand(2,8)-0.5
theta2=np.random.rand(8,1)-0.5

# Necessary - the bias terms should have same number of dimensions
# as the layer. For some reason you had one bias vector per example.
# (You could still use np.random.rand(8) and np.random.rand(1))
b1=np.zeros(8)
b2=np.zeros(1)

alpha=0.01
# Regularisation not necessary for XOR, because you have a complete training set.
# You could have lamda=0.0, but I have left a value here just to show it works.
lamda=0.001


Training - Forward propagation

# More iterations than you might think! This is because we have
# so little training data, we need to repeat it a lot.
for i in range(1,40000):
    z1=x.dot(theta1)+b1
    h1=1/(1+np.exp(-z1))
    z2=h1.dot(theta2)+b2
    h2=1/(1+np.exp(-z2))


Training - Back propagation

    # This dz term assumes binary cross-entropy loss
    dz2 = h2-y 
    # You could also have stuck with squared error loss, the extra h2 terms
    # are the derivative of the sigmoid transfer function. 
    # It converges slower though:
    # dz2 = (h2-y) * h2 * (1-h2)

    # This is just the same as you had before, but with less temp variables
    dw2 = np.dot(h1.T, dz2)
    db2 = np.sum(dz2, axis=0)

    # The derivative of sigmoid is h1 * (1-h1), NOT dh1*(1-dh1)
    dz1 = np.dot(dz2, theta2.T) * h1 * (1-h1)
    dw1 = np.dot(x.T, dz1)
    db1 = np.sum(dz1, axis=0)

    # The L2 regularisation terms ADD to the gradients of the weights
    dw2 += lamda * theta2
    dw1 += lamda * theta1

    theta1 += -alpha * dw1
    theta2 += -alpha * dw2

    b1 += -alpha * db1
    b2 += -alpha * db2


Prediction

This is where you can kick yourself, you forgot to use the biases!

input1=np.array([[0,0],[1,1],[0,1],[1,0]])
z1=np.dot(input1,theta1)+b1
h1=1/(1+np.exp(-z1))
z2=np.dot(h1,theta2)+b2
h2=1/(1+np.exp(-z2))

print(h2)


When I run the above code I get a correct-looking output

[[ 0.01031446]
 [ 0.0201576 ]
 [ 0.9824826 ]
 [ 0.98584079]]




In summary your three big errors were the wrong dimension for the bias vectors in setup, incorrect derivatives for sigmoid function (using correct form, but with wrong variable) and forgetting to use bias at all when predicting at the end. Other details are still worth noting, but would not have prevented you getting something working.