Let's assume that you are given a model simulating the number of events as a function of a certain binned variable. The simulation also depends on a parameter that is normally distributed with zero mean and unit sigma.
$$N_i(X) := \text{number of events in the }i\text{-bin as a function of the }X\text{ parameter normally distributed with }\mu=0\text{ and }\sigma=1$$
By running the simulation many times we can measure the average dependence of the number of events on $X$, where in the simplest case the relationship is approximately linear $N_i(X) \approx b_i \cdot X$. In such a case we call the $X$ parameter "gaussian parameter".
If the relationship is not linear, we call the parameter "non-gaussian" and we sample the average number of events at $X=0,\pm 1, \pm 2, ... \pm n$ and approximate $N_i()$ as a 3-order spline.
My problem is that I have to calculate the covariance between the $i$-bin and $j$-bin in case the parameter is "non-gaussian" but I am stuck.
I can calculate the covariance when the parameter is gaussian and the result is:
$$1)~~Cov_{ij}[X]:=E[(N_i(X)-E[N_i(X)])(N_j(X)-E[N_j(X)])]=E[N_i(X)N_i(X)]-E[N_i(X)]E[N_j(X)]=$$ $$= E[b_i b_j X^2] - \bar{N}_i\bar{N}_j=b_i b_jE[X^2] - \bar{N}_i\bar{N}_j=N_i^{+}N_j^{+} - \bar{N}_i\bar{N}_j$$
where I used the fact that $E[X^2]=\sigma^2$=1 and I defined $N_i^{+}:=N_i(X=+1)$ and $\bar{N}_i := N_i(X=0)$. I use this notation because it is conventional in our group and because it is easier to generalize. I cross-checked this equation with other people and it should be correct.
I would be very grateful if you could point me in the right direction of how to generalize equation 1) in case the relationship between $N_i$ and $X$ is generic and all I have are the sampled points $N_i(X=0,\pm 1,...)$.
Thank you.