Friends, I want to share the details of my small algorithm to correct the moments of arrays of standard normals associated with equiprobable normal grid subdivisions. These arrays of standard normals would later be used for polynomial regression on the data to find parameters of Z-series from the data.

Suppose there are N data members and we want to do a hermite polynomial regression on this data to find Z-series parameters of data random variable. We will divide the standard normal density into N subdivisions with each subdivision sharing an equal probability mass of 1/N.

However discrete values of Zn when integrated on probability mass 1/N, can have moments different from moments of standard normal density which will create inaccuracies in our polynomial regression to find Z-series parameters of the data random variable.

I thought of multiplying the data array with an expression comprising sum of hermite polynomials with unspecified coefficients. We would write equations of moments and then find coefficients of hermite polynomials that pin down the moments of new data precisely equal to target moments.

One consideration is that our Zn array is perfectly symmetric around zero so all odd moments are zero. We want to add only even hermite polynomials in the expression of hermite polynomials with unspecified coefficients(that is multiplied to original Zn array to make its moments equal to target moments). We do not want odd hermite polynomials in this expression with unspecified coefficients as they would introduce skew in our array of equi-probable standard normal values(and then we would have to do unnecessary work to cancel the skew). In order to have zero odd moments(as in case of standard normals), we only multiply with a series in even hermite polynomials with unspecified coefficients that preserves perfect symmetry around zero so we would not have to bother with any odd moments and all odd moments will be zero by construction due to symmetry of our Zn array around zero.

Our original standard normal Zn array associated with equiprobable grid is multiplied with a series in hermite polynomials as

[$]Z_{n,Corrected} \, = \, Z_n \, (\alpha \, + \beta \, (Z_n^2-1) \, + \, \gamma \, ({Z_n}^4 \, + \, 6 \, {Z_n}^2 \, + \, 3 \,))\,[$]

We find coefficients [$]\alpha[$], [$]\beta[$], and [$]\gamma[$] so that our objective function matching moments of [$]Z_{n,Corrected} \, [$] with target moments is maximized.

[$]Moment2\, = \sum\limits_{n=1}^{n=N} \, {Z_{n,Corrected}}^2 \, \frac{1}{N} [$]

[$]= \, \sum\limits_{n=1}^{n=N} \, {Z_n}^2 \, {(\alpha \, + \beta \, (Z_n^2-1) \, + \, \gamma \, ({Z_n}^4 \, + \, 6 \, {Z_n}^2 \, + \, 3 \,))}^2\, \frac{1}{N} [$]

[$]Moment4\, = \sum\limits_{n=1}^{n=N} \, {Z_{n,Corrected}}^4 \, \frac{1}{N} [$]

[$]= \, \sum\limits_{n=1}^{n=N} \, {Z_n}^4 \, {(\alpha \, + \beta \, (Z_n^2-1) \, + \, \gamma \, ({Z_n}^4 \, + \, 6 \, {Z_n}^2 \, + \, 3 \,))}^4\, \frac{1}{N} [$]

[$]Moment6\, = \sum\limits_{n=1}^{n=N} \, {Z_{n,Corrected}}^6 \, \frac{1}{N} [$]

[$]= \, \sum\limits_{n=1}^{n=N} \, {Z_n}^6 \, {(\alpha \, + \beta \, (Z_n^2-1) \, + \, \gamma \, ({Z_n}^4 \, + \, 6 \, {Z_n}^2 \, + \, 3 \,))}^6\, \frac{1}{N} [$]

Here is an iterative optimization matlab program that finds [$]Z_{n,Corrected}[$] whose moments match the standard normal moments with great precision. I improved this program and now it takes only 3-7 seconds as opposed to 30-40 secs it was taking yesterday. So if great speed is not the concern, you can easily embed it in your programs.

.

.

.

```
function [ZnNew] = MatchMomentsOfIdealZ(Zn)
Ndata=length(Zn);
if(Ndata<=100)
Iterations=3500;
end
if((Ndata>100)&&(Ndata<=500))
Iterations=2000;
end
if((Ndata>500)&&(Ndata<=1000))
Iterations=1500;
end
if((Ndata>1000)&&(Ndata<=3000))
Iterations=1000;
end
if((Ndata>3000)&&(Ndata<=6000))
Iterations=750;
end
if((Ndata>6000)&&(Ndata<=10000))
Iterations=550;
end
if(Ndata>10000)
Iterations=50;
end
alpha=1;
beta=0;
gamma=0;
ZnNew(1:Ndata)=Zn(1:Ndata).*(alpha+beta*(Zn(1:Ndata).^2-1)+gamma*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata
Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata
Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata
w1=5000;
w2=200;
ObjFuncBest=(Moment2-1).^2*w1+(Moment4-3).^2*w2+(Moment6-15).^2;
delta_alpha=.00025;
delta_beta=.0002;
delta_gamma=.0002;
alphaSuccessCount=0;
betaSuccessCount=0;
gammaSuccessCount=0;
SuccessPrevalpha=0;
SuccessPrevbeta=0;
SuccessPrevgamma=0;
mul1=.9;
mul2=1.1;
for iter=1:Iterations
alphaNew=alpha+delta_alpha;
betaNew=beta;
gammaNew=gamma;
ZnNew(1:Ndata)=Zn(1:Ndata).*(alphaNew+betaNew*(Zn(1:Ndata).^2-1)+gammaNew*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata;
Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata;
Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata;
ObjFuncNew=(Moment2-1).^2*w1+(Moment4-3).^2*w2+(Moment6-15).^2;
if(ObjFuncNew<ObjFuncBest)
alpha=alphaNew;
ObjFuncBest=ObjFuncNew;
if(SuccessPrevalpha==1)
alphaSuccessCount=alphaSuccessCount+1;
end
SuccessPrevalpha=1;
else
alphaNew=alpha-delta_alpha;
betaNew=beta;
gammaNew=gamma;
ZnNew(1:Ndata)=Zn(1:Ndata).*(alphaNew+betaNew*(Zn(1:Ndata).^2-1)+gammaNew*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata;
Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata;
Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata;
ObjFuncNew=(Moment2-1).^2*w1+(Moment4-3).^2*w2+(Moment6-15).^2;
if(ObjFuncNew<ObjFuncBest)
alpha=alphaNew;
ObjFuncBest=ObjFuncNew;
if(SuccessPrevalpha==1)
alphaSuccessCount=alphaSuccessCount+1;
end
SuccessPrevalpha=1;
else
SuccessPrevalpha=0;
alphaSuccessCount=0;
delta_alpha=delta_alpha*mul1;
end
end
alphaNew=alpha;
betaNew=beta+delta_beta;
gammaNew=gamma;
ZnNew(1:Ndata)=Zn(1:Ndata).*(alphaNew+betaNew*(Zn(1:Ndata).^2-1)+gammaNew*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata;
Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata;
Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata;
ObjFuncNew=(Moment2-1).^2*w1+(Moment4-3).^2*w2+(Moment6-15).^2;
if(ObjFuncNew<ObjFuncBest)
beta=betaNew;
ObjFuncBest=ObjFuncNew;
if(SuccessPrevbeta==1)
betaSuccessCount=betaSuccessCount+1;
end
SuccessPrevbeta=1;
else
alphaNew=alpha;
betaNew=beta-delta_beta;
gammaNew=gamma;
ZnNew(1:Ndata)=Zn(1:Ndata).*(alphaNew+betaNew*(Zn(1:Ndata).^2-1)+gammaNew*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata;
Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata;
Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata;
ObjFuncNew=(Moment2-1).^2*w1+(Moment4-3).^2*w2+(Moment6-15).^2;
if(ObjFuncNew<ObjFuncBest)
beta=betaNew;
ObjFuncBest=ObjFuncNew;
if(SuccessPrevbeta==1)
betaSuccessCount=betaSuccessCount+1;
end
SuccessPrevbeta=1;
else
SuccessPrevbeta=0;
betaSuccessCount=0;
delta_beta=delta_beta*mul1;
end
end
alphaNew=alpha;
betaNew=beta;
gammaNew=gamma+delta_gamma;
ZnNew(1:Ndata)=Zn(1:Ndata).*(alphaNew+betaNew*(Zn(1:Ndata).^2-1)+gammaNew*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata;
Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata;
Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata;
ObjFuncNew=(Moment2-1).^2*w1+(Moment4-3).^2*w2+(Moment6-15).^2;
if(ObjFuncNew<ObjFuncBest)
gamma=gammaNew;
ObjFuncBest=ObjFuncNew;
if(SuccessPrevgamma==1)
gammaSuccessCount=gammaSuccessCount+1;
end
SuccessPrevgamma=1;
else
alphaNew=alpha;
betaNew=beta;
gammaNew=gamma-delta_gamma;
ZnNew(1:Ndata)=Zn(1:Ndata).*(alphaNew+betaNew*(Zn(1:Ndata).^2-1)+gammaNew*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata;
Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata;
Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata;
ObjFuncNew=(Moment2-1).^2*w1+(Moment4-3).^2*w2+(Moment6-15).^2;
if(ObjFuncNew<ObjFuncBest)
gamma=gammaNew;
ObjFuncBest=ObjFuncNew;
if(SuccessPrevgamma==1)
gammaSuccessCount=gammaSuccessCount+1;
end
SuccessPrevgamma=1;
else
delta_gamma=delta_gamma*mul1;
SuccessPrevgamma=0;
gammaSuccessCount=0;
end
end
if(alphaSuccessCount>4)
delta_alpha=delta_alpha*mul2;
end
if(betaSuccessCount>4)
delta_beta=delta_beta*mul2;
end
if(gammaSuccessCount>4)
delta_gamma=delta_gamma*mul2;
end
end
ZnNew(1:Ndata)=Zn(1:Ndata).*(alpha+beta*(Zn(1:Ndata).^2-1)+gamma*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata;
Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata;
Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata;
Moment2
Moment4
Moment6
%str=input("Look at moment values")
end
```