Here is my code, Nintegrate of mixture gaussian distribution:
Gaussian[x_?NumericQ, y_?NumericQ, mu1_?NumericQ, mu2_?NumericQ, sigma_?NumericQ]:= PDF[MultinormalDistribution[{mu1, mu2}, {{sigma, 0}, {0, sigma}}], {x, y}]
q[x_?NumericQ, y_?NumericQ] := 0.1 * Gaussian[x,y, -0.5, 0.5, 1] + 0.2 * Gaussian[x,y, 0.5, 0.5, 1] + 0.3 * Gaussian[x,y, 0.5, -0.5, 1] +
0.4 * Gaussian[x,y, -0.5, -0.5, 1]
p[x_?NumericQ, y_?NumericQ, w1_?NumericQ, w2_?NumericQ, w3_?NumericQ, mu11_?NumericQ, mu12_?NumericQ, mu21_?NumericQ, mu22_?NumericQ, mu31_?NumericQ, mu32_?NumericQ, mu41_?NumericQ, mu42_?NumericQ, sigma_?NumericQ] := w1 * Gaussian[x, y, mu11, mu12, sigma] + w2 * Gaussian[x,y, mu21, mu22, sigma] + w3 * Gaussian[x,y, mu31, mu32, sigma] +
(1 - w1 - w2 - w3) * Gaussian[x,y, mu41, mu42, sigma]
K[w1_?NumericQ, w2_?NumericQ, w3_?NumericQ, mu11_?NumericQ, mu12_?NumericQ, mu21_?NumericQ, mu22_?NumericQ, mu31_?NumericQ, mu32_?NumericQ, mu41_?NumericQ, mu42_?NumericQ, sigma_?NumericQ] := NIntegrate[q[x, y] * Log[q[x, y]/p[x, y, w1, w2, w3, mu11, mu12, mu21, mu22, mu31, mu32, mu41, mu42, sigma]], {x, -Infinity, Infinity}, {y, -Infinity, Infinity}]
K[0.1,0.2,0.2,1,1, 1,-1, -1,1, -1,-1, 2]
Kesi[z_?NumericQ] = NIntegrate[(K[w1, w2, w3, mu11, mu12, mu21, mu22, mu31, mu32, mu41, mu42, sigma]) ^ z, {mu11, -1, 1}, {mu12, -1, 1}, {mu21, -1, 1}, {mu22, -1, 1}, {mu31, -1, 1}, {mu32, -1, 1}, {mu41, -1, 1}, {mu42, -1, 1}, {sigma, 0.5, 1.5}, {w1, w2, w3} \[Element] RegionIntersection[ImplicitRegion[w1+w2+w3 <= 1, {w1, w2, w3}], ImplicitRegion[w1>= 0, {w1, w2, w3}], ImplicitRegion[w2>=0, {w1, w2,w3}], ImplicitRegion[w3 >=0, {w1, w2, w3}]]]
Kesi[1]