#include <iostream>
#include <fstream>
#include <Eigen/Dense>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
using
namespace
Eigen;
using
namespace
std;
#define PI 3.141592654
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
double
ran1(
long
*idum)
{
int
j;
long
k;
static
long
iy=0;
static
long
iv[NTAB];
double
temp;
if
(*idum <= 0 || !iy) {
if
(-(*idum) < 1) *idum=1;
else
*idum = -(*idum);
for
(j=NTAB+7;j>=0;j--) {
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if
(*idum < 0) *idum += IM;
if
(j < NTAB) iv[j] = *idum;
}
iy=iv[0];
}
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if
(*idum < 0) *idum += IM;
j=iy/NDIV;
iy=iv[j];
iv[j] = *idum;
if
((temp=AM*iy) > RNMX)
return
RNMX;
else
return
temp;
}
double
gasdev(
long
*idum)
{
double
ran1(
long
*idum);
static
int
iset=0;
static
double
gset;
double
fac,rsq,v1,v2;
if
(*idum<0) iset=0;
if
(iset == 0){
do
{
v1=2.0*ran1(idum)-1.0;
v2=2.0*ran1(idum)-1.0;
rsq=v1*v1+v2*v2;
}
while
(rsq >= 1.0 || rsq == 0.0);
fac=
sqrt
(-2.0*
log
(rsq)/rsq);
gset=v1*fac;
iset=1;
return
v2*fac;
}
else
{
iset=0;
return
gset;
}
}
double
gamdev(
int
ia,
long
*idum)
{
double
ran1(
long
*idum);
int
j;
double
am,e,s,v1,v2,x,y;
if
(ia<6) {
x=1.0;
for
(j=1;j<=ia;j++) x*=ran1(idum);
x=-
log
(x);
}
else
{
do
{
do
{
do
{
v1=ran1(idum);
v2=2.0*ran1(idum)-1.0;
}
while
(v1*v1+v2*v2>1.0);
y=v2/v1;
am=ia-1;
s=
sqrt
(2.0*am+1.0);
x=s*y+am;
}
while
(x<=0.0);
e=(1.0+y*y)*
exp
(am*
log
(x/am)-s*y);
}
while
(ran1(idum)>e);
}
return
x;
}
double
P(
int
i,
int
j)
{
double
ans=1;
int
i1;
for
(i1=1;i1<=j;i1++)
ans=ans*i1;
for
(i1=1;i1<=j-i;i1++)
ans=ans/
double
(i1);
return
ans;
}
int
main()
{
int
n,p,i,j,l,rep,v,k0,q,qnum;
n=50;
p=50;
rep=1000;
k0=5;
qnum=49;
MatrixXd G,z(n,p),x(n,p),sum(rep,(qnum+1)),sigma(rep,(qnum+1)),O,C,xcolsum,xmean;
G = MatrixXd::Identity(p,p);
O = MatrixXd::Ones(n,n);
for
(j=1;j<(k0+1);j++)
{
for
(i=0;i<p-j;i++)
{
G(i+j,i) = 0.4;
}
}
clock_t
tStart =
clock
();
for
(v=0;v<rep;v++)
{
long
idum=-(v+12345);
for
(i=0;i<n;i++)
{
for
(j=0;j<p;j++)
{
z(i,j)=gasdev(&idum);
}
}
x=z*G.transpose();
xcolsum = x.colwise().sum();
xmean = O*x/n;
C=(x-O*x/n).transpose()*(x-O*x/n)/(n-1);
for
(q=0;q<(qnum+1);q++)
{
MatrixXd Y(n,(p-q));
for
(i=0;i<n;i++)
{
for
(l=0;l<(p-q);l++)
{
Y(i,l) = (x(i,l)-xmean(l))*(x(i,(l+q))-xmean(l+q)) - C(l,(l+q));
}
}
MatrixXd YY(n,n),YYcolsum;
YY=Y*Y.transpose();
double
sigma1,sum1,temp1,temp2,temp3;
sigma1 = 0;
sum1=0;
temp1=0;
temp2=0;
temp3=0;
for
(i=0;i<n;i++)
{
for
(j=0;j<i;j++)
{
for
(l=0;l<p-q;l++)
{
temp1 += x(i,l)*x(i,l+q)*x(j,l)*x(j,l+q);
temp2 += (x(i,l)*x(j,l)*x(j,l+q)*x(j,l+q) - x(i,l)*x(j,l)*x(j,l+q)*xcolsum(l+q));
temp3 += (x(i,l)*x(j,l+q)*xcolsum(l)*xcolsum(l+q) + x(i,l)*x(i,l)*x(j,l+q)*x(j,l+q) - x(i,l)*x(i,l)*x(j,l+q)*xcolsum(l+q) - x(i,l)*x(j,l+q)*x(j,l+q)*xcolsum(l));
}
}
for
(j=(i+1);j<n;j++)
{
for
(l=0;l<p-q;l++)
{
temp1 += x(i,l)*x(i,l+q)*x(j,l)*x(j,l+q);
temp2 += (x(i,l)*x(j,l)*x(j,l+q)*x(j,l+q) - x(i,l)*x(j,l)*x(j,l+q)*xcolsum(l+q));
temp3 += (x(i,l)*x(j,l+q)*xcolsum(l)*xcolsum(l+q) + x(i,l)*x(i,l)*x(j,l+q)*x(j,l+q) - x(i,l)*x(i,l)*x(j,l+q)*xcolsum(l+q) - x(i,l)*x(j,l+q)*x(j,l+q)*xcolsum(l));
}
}
sigma1 += YY(i,i)*YY(i,i);
}
sum(v,q)=temp1/n/(n-3)+temp2*(2*n-3)/n/(n-1)/(n-2)/(n-3)+temp3/n/(n-1)/(n-2)/(n-3);
sigma(v,q)=((YY*YY).trace()-sigma1)/n/(n-1);
}
}
ofstream f;
f.open (
"Sigma_n50p50normal.txt"
);
f << sigma;
f.close();
ofstream f1;
f1.open (
"Dnq_n50p50normal.txt"
);
f1 << sum;
f1.close();
return
0;
}