Hello Frank,
Thank you very much for your response. I was able to write a function as follows:
from TCN.poly_music.myintegrate2 import integrand, integrand_alpha, integrand_beta, integrand_mu, integrand_sigma
from scipy.integrate import quad 
import torch

class my_loss(torch.autograd.Function):
    @staticmethod
    def forward(ctx,a,b,m,s,a2,b2,m2,s2):
        result = quad(integrand, 0, 5*10**-7, args=(a,b,m,s,a2,b2,m2,s2), points=(max(0,min(m-2*s,m2-2*s2)), max(m+2*s,m2+2*s2)))[0]
        ctx.save_for_backward(a,b,m,s,a2,b2,m2,s2)
        return torch.FloatTensor([result])


    @staticmethod
    def backward(ctx, grad_output):
        a,b,m,s,a2,b2,m2,s2 = ctx.saved_variables
        grad_alpha = quad(integrand_alpha, 0, 5*10**-7, args=(a,b,m,s,a2,b2,m2,s2), points=(max(0,min(m-2*s,m2-2*s2)), max(m+2*s,m2+2*s2)))[0]
        grad_beta = quad(integrand_beta, 0, 5*10**-7, args=(a,b,m,s,a2,b2,m2,s2), points=(max(0,min(m-2*s,m2-2*s2)), max(m+2*s,m2+2*s2)))[0]
        grad_mu = quad(integrand_mu, 0, 5*10**-7, args=(a,b,m,s,a2,b2,m2,s2), points=(max(0,min(m-2*s,m2-2*s2)), max(m+2*s,m2+2*s2)))[0]
        grad_sig = quad(integrand_sigma, 0, 5*10**-7, args=(a,b,m,s,a2,b2,m2,s2), points=(max(0,min(m-2*s,m2-2*s2)), max(m+2*s,m2+2*s2)))[0]
            
    
        return torch.FloatTensor(grad_alpha), torch.FloatTensor(grad_beta), torch.FloatTensor(grad_mu), torch.FloatTensor(grad_sig), None

In this function integrand, integrand_alpha, integrand_beta, integrand_mu and integrand_sigme are function of this sort:
def integrand(y,a,b,m,s,a2,b2,m2,s2): 

     return abs(((1/2)*b*math.exp(1)**((-1....)

And the  first values that I get are as follows:
     grad_alpha = 0.210814
     grad_beta = -0.000110
     grad_mu =   563428.64
     grad_sig =   3272315.65

In the main program I used this one:
    loss = my_loss.apply(a[0],b[0],m[0],s[0],a2,b2,m2,s2)

Here the loss is of this sort:
    tensor([1.1734])

The rest of the code is the same as I had written in the original question. I got an error as this one:
    element 0 of tensors does not require grad and does not have a grad_fn

Later I checked in Internet and I saw that this error was related to no gradients in the tensor and a quick fix would be to use Variable after the definition of loss.
        loss = my_loss.apply(a[0],b[0],m[0],s[0],a2,b2,m2,s2)
        loss = Variable(loss, requires_grad = True)

This lets the code run but the results show that the loss function does not decrease
 Epoch  1 | lr 0.00100 | loss 1.62780
 Epoch  1 | lr 0.00100 | loss 1.63068
 Epoch  1 | lr 0.00100 | loss 1.61624
 Epoch  1 | lr 0.00100 | loss 1.62332
 Epoch  1 | lr 0.00100 | loss 1.62789
 Epoch  1 | lr 0.00100 | loss 1.64436

I think there is something missing about the Function that I defined, else I should have seen the loss function decreasing. To be on the safe side, I also checked the code for a single parameter and without custom loss function, i.e. with MSE loss that I defined over torch. This worked very well and gave some expected results. I calculated the gradients with Mathematica and tested them in Python and got the same results with Matlab. So the gradients should be OK.
What should then be missing with my Function or the appliaction of it in the main code?