Question: How to obtain an analytical solution with respect to the parameters of a parametric equation

Questions are as follows:

Consider a g function, for example: g = 1 - (x_1 - 1)^2/9 - (x_2 - 1)^3/16;

Among them, x_1 and x_2 each conform to the normal distribution, and each has a mean mu_1 mu_2 and a variance sigma_1 sigma_2;

mu_1 and mu_2 both obey uniform distribution (-0.2, 0.2) respectively; sigma_1 and sigma_2 both obey uniform distribution (0.8, 1.2) respectively;

Set parameter theta = (theta_1, theta_2, theta_3, theta_4)^T = (mu_1, mu_2, sigma_1, sigma_2)^T;

It should be possible to obtain an analytical solution with only theta for the mean and variance of the g function, such as:

1. Ey = -(1/9)*theta[1]^2 + (2/9)*theta[1] + 137/144 - (1/16)*theta[2]^3 + (3/16)* theta[2]^2 - (3/16)*theta[2] - (3/16*(-1+theta[2]))*theta[4]^2 - (1/9)*theta[3 ]^2;
2. Vy = (1/20736)*(9*theta[2]^3 + 16*theta[1]^2 - 27*theta[2]^2 - 32*theta[1] + 27*theta[2] - 137)^2 + (15/256)*theta[4]^6 + (45/256*(theta[2]^2 - 2*theta[2] + 1))*theta[4]^4 + ( 1/768*(45*theta[2]^4 + 32*theta[1]^2*theta[2] - 180*theta[2]^3 - 32*theta[1]^2 - 64*theta[ 1]*theta[2] + 270*theta[2]^2 + 64*theta[1] - 436*theta[2] + 301))*theta[4]^2 + (1/27)*theta[ 3]^4 + (1/216)*theta[3]^2*(3*theta[2]^3 + 16*theta[1]^2 - 9*theta[2]^2 - 32*theta[ 1] + 9*theta[2] - 35) + (1/24)*theta[3]^2*(-1+theta[2])*theta[4]^2;

But I can't find the above analytical solution. The code I wrote in maple is as follows:

restart;

with(Statistics);

x[1] := RandomVariable(Normal(theta[1], theta[3]));

x[2] := RandomVariable(Normal(theta[2], theta[4]));

g := 1 - (x[1] - 1)^2/9 - (x[2] - 1)^3/16;

theta[1] := RandomVariable(Uniform(-0.2, 0.2));

theta[2] := RandomVariable(Uniform(-0.2, 0.2));

theta[3] := RandomVariable(Uniform(0.8, 1.2));

theta[4] := RandomVariable(Uniform(0.8, 1.2));

mean := Mean(g);

variance := Variance(g);

mean;

variance;

As a result, after clicking Run, the following problems occur:

Problem 1: The program takes a very long time to run

Problem 2: The program results did not meet expectations, because the results I expected were：

Ey = -(1/9)*theta[1]^2 + (2/9)*theta[1] + 137/144 - (1/16)*theta[2]^3 + (3/16)* theta[2]^2 - (3/16)*theta[2] - (3/16*(-1+theta[2]))*theta[4]^2 - (1/9)*theta[3 ]^2;

Vy = (1/20736)*(9*theta[2]^3 + 16*theta[1]^2 - 27*theta[2]^2 - 32*theta[1] + 27*theta[2] - 137)^2 + (15/256)*theta[4]^6 + (45/256*(theta[2]^2 - 2*theta[2] + 1))*theta[4]^4 + ( 1/768*(45*theta[2]^4 + 32*theta[1]^2*theta[2] - 180*theta[2]^3 - 32*theta[1]^2 - 64*theta[ 1]*theta[2] + 270*theta[2]^2 + 64*theta[1] - 436*theta[2] + 301))*theta[4]^2 + (1/27)*theta[ 3]^4 + (1/216)*theta[3]^2*(3*theta[2]^3 + 16*theta[1]^2 - 9*theta[2]^2 - 32*theta[ 1] + 9*theta[2] - 35) + (1/24)*theta[3]^2*(-1+theta[2])*theta[4]^2; 

so,i dont know where i am wrong ,how could i obtain the expected solution?