How to integrate directional data

I have a data vector that consists of directions measured in radians. I wish to create a probability density function. As the data is circular and multimodal I use a kernel density estimate with von Mises distribution as the basis function. I fit a von Mises function to each data point and sum the results to obtain a smooth distribution. To obtain a probability density I simply divide each point in the distribution by the integral of the whole distribution. However, my results seem odd after the integration as the maximum value in the pdf is larger than 1. I think it might be related to how I do the integration, I use the trapz command but I am not sure if this is appropriate for circular data.