C++ Program

#include<iostream>

#include<fstream>

#include<cmath>

#include<list>

#include "Unode.h"

#include "Vnode.h"

using namespace std;

void computePVal(double* t, double* count, int numt, double tobs);

void getObservables(int I, int* J, int** n, int** z, double** x, double& tobs, int&

s1obs, int& s2obs, int& s3obs);

//Functions for U

void insertU(list<Unode>* theList, Unode& entry);

void printU(list<Unode>& theList);

void minFeasibleValsU(int I, int* J, int** n, int inow, int jnow, int s1obs, int

s2obs, int& s1min, int& s2min);

//Functions for V should mirror those for U

void insertV(list<Vnode>* theList, Vnode& entry);

void printU(list<Vnode>& theList);

void minFeasibleValsV(int I, int* J, int** n, int inow, int s1obs,

int s2obs, int s3obs, int& s1min, int& s2min, int& s3min);

//General function but needs a rewrite

void computeDistribution(int I, int* J, int** n, double** x, int s1obs, int s2obs,

int s3obs, double tobs);

int main(void)

{

//observed values of the sufficient statistics

double tobs = -1;

int s1obs = -1;

int s2obs = -1;

int s3obs = -1;

//the total number of primary clusters

int I = 1;

cout << "Enter the number of primary clusters: ";

cin >> I;

//arrays for the various data

int* J = new int[I]; //Vector with the number of secondary clusters in each

primary cluster

int** n = new int*[I]; //number of observations per cluster

int** z = new int*[I]; //number of successes observed per cluster

double** x = new double*[I]; //value of the covariate in each secondary cluster

//To make the algorithm efficient, these should

take values from

//a relatively small discrete set (integers,

half-integers, quarter-ints, etc)

getObservables(I, J, n, z, x, tobs, s1obs, s2obs, s3obs);

computeDistribution(I, J, n, x, s1obs, s2obs, s3obs, tobs);

return 0;

}//end main

void getObservables(int I, int* J, int** n, int** z, double** x, double& tobs, int&

s1obs, int& s2obs, int& s3obs)

{

for(int i=0; i<I; i++)

{

cout<< "How many secondary clusters are there in primary cluster " << i <<

"?"<< endl;

cin >> J[i];

n[i]=new int[J[i]];

z[i]=new int[J[i]];

x[i]=new double[J[i]];

for(int j=0; j<J[i]; j++)

{

cout << "In secondary cluster " << i << "," << j << ", enter the following"

<< endl;

cout << "n: ";

cin >> n[i][j];

cout << "z: ";

cin >> z[i][j];

cout << "x: ";

cin >> x[i][j];

}

}

tobs = 0;

s1obs = 0;

s2obs = 0;

s3obs = 0;

for(int i=0; i<I; i++)

{

int nTot = 0;

int zTot = 0;

for(int j=0; j<J[i]; j++)

{

s1obs += z[i][j];

s2obs += z[i][j]*(n[i][j]-z[i][j]);

tobs += z[i][j]*x[i][j];

nTot += n[i][j];

zTot += z[i][j];

}

s3obs += zTot*(nTot-zTot);

}

}

void computePVal(double* t, double* count, int numt, double tobs)

{

double total = 0, sum=0;

for(int ndx = 0; ndx<numt; ndx++)

{

total += count[ndx];

if(t[ndx]>=tobs)

sum+=count[ndx];

}

double pval = sum/total;

cout << "One-sided Pvalue is: " << pval << endl;

}

//***************FUNCTIONS INVOLVING JUST U.

void insertU(list<Unode>* theList, Unode& entry)

{

if((*theList).begin()==(*theList).end()) //A blank list. This case is not

tested in the demo.

{

Unode toInsert(entry);

(*theList).push_back(toInsert);

}

else

{

list<Unode>::iterator iterU;

iterU = (*theList).begin();

while(entry > *iterU && iterU != (*theList).end())

    iterU++;

if(entry == *iterU)

iterU->augmentCount(entry);

else //entry belongs in front of *iterU...this works even

if iterU is the end.

{

Unode toInsert(entry);

(*theList).insert(iterU, toInsert);

}

}

}

{

    for(list<Unode>::iterator iterU=theList.begin(); iterU != theList.end();

++iterU)

        (*iterU).print();

}

void minFeasibleValsU(int I, int* J, int** n, int inow, int jnow, int s1obs,

int s2obs, int& s1min, int& s2min)

{

int s1poss = 0, s2poss = 0;

//sum the possible values after the current one in this cluster

for(int j = jnow+1; j<J[inow]; j++)

{

s1poss += n[inow][j]; //the largest possible addition to s1 in that cluster

s2poss += (n[inow][j]*n[inow][j])/4; //the largest possible addition to s2

}

//sum the possible values in later clusters

for(int i=inow+1; i<I; i++)

{

for(int j=0; j<J[i]; j++)

{

s1poss += n[i][j]; //the largest possible addition to s1 in that cluster

s2poss += (n[i][j]*n[i][j])/4; //the largest possible addition to s2

}

}

//these are the minimum values at this stage that would make it possible

to eventually reach observeds.

s1min = s1obs - s1poss;

s2min = s2obs - s2poss;

}

//***************FUNCTIONS INVOLVING JUST V.

void insertV(list<Vnode>* theList, Vnode& entry)

{

if((*theList).begin()==(*theList).end())

{

Vnode toInsert(entry);

(*theList).push_back(toInsert);

}

else

{

list<Vnode>::iterator iterV;

iterV = (*theList).begin();

while(entry > *iterV && iterV != (*theList).end())

    iterV++;

if(entry == *iterV)

(*iterV).augmentCount(entry);

else         //entry belongs in front of *iterU...this works even if

iterV is the end.

{

Vnode toInsert(entry);

(*theList).insert(iterV, toInsert);

}

}

}//end insertV

void printV(list<Vnode>& theList)

{

    for(list<Vnode>::iterator iterU=theList.begin(); iterU != theList.end();

++iterU)

(*iterU).print();

}

void minFeasibleValsV(int I, int* J, int** n, int inow, int s1obs,

    int s2obs, int s3obs, int& s1min, int& s2min, int& s3min)

{

//note: we don't need a preliminary loop this time since we are at the end

of the primary cluster

int s1poss = 0, s2poss = 0, s3poss = 0;

//sum the possible values in later primary clusters

for(int i=inow+1; i<I; i++)

{

//Figure out the largest possible addition to s3 from primary cluster i.

int nTot = 0;

for(int j=0; j<J[i]; j++)

nTot += n[i][j];

s3poss += (nTot*nTot)/4;

//Sum the possible values for secondary clusters within the primary cluster

for(int j=0; j<J[i]; j++)

{

s1poss += n[i][j];

s2poss += (n[i][j]*n[i][j])/4;

}

}

//these are the minimum values at this stage that would make it

possible to eventually reach observeds.

s1min = s1obs - s1poss;

s2min = s2obs - s2poss;

s3min = s3obs - s3poss;

}//end getObservables

//***********************Compute Distribution.

void computeDistribution(int I, int* J, int** n, double** x, int s1obs, int

s2obs, int s3obs, double tobs)

{

//defining V lists and iterator for scanning

list<Vnode> newV;

list<Vnode> oldV;

list<Vnode>::iterator iterV;

//making a base V list before i loop.

Vnode addV;

insertV(&newV, addV);

for(int i=0; i< I; i++)

{

//making the copy of the V list for next iteration

oldV = newV;

newV = list<Vnode>();

//defining U lists and iterator for scanning

list<Unode> newU;

list<Unode> oldU;

list<Unode>::iterator iterU;

//making a base U list before j loop.

Unode addU;

insertU(&newU, addU);

for(int j=0; j<J[i]; j++)

{

//copying the U list for next iteration

oldU = newU;

newU = list<Unode>();

for(int z=0; z<=n[i][j]; z++)

{

//Creating the set of additions to suff statistics if

z[i][j]=z were true

addU = Unode(n[i][j],z,x[i][j]);

//Adding those additions to existing suff statistics

for(iterU = oldU.begin(); iterU != oldU.end(); iterU++)

{

//the new value if z[i][j]=z after adding

Unode addedU = addU + (*iterU);

//feasibility check...can't easily do minimums here,

so only check maximums

if(addedU.isFeasible(s1obs, s2obs))

insertU(&newU, addedU);

}

}//end z-loop

//releasing memory, since we have now finished with the i,j-th

second-stage cluster

oldU.clear();

}//End j loop...Unew is now complete for the ith primary cluster

/****Updating newV using oldV and newV****/

//In order to update, we need the total number of observations in

i-th primary cluster

int nTot = 0;

for(int j=0; j<J[i]; j++)

nTot += n[i][j];

//Loop, one iteration for each element of newU

for(iterU = newU.begin(); iterU != newU.end(); iterU++)

{

//create that element

Vnode Vadd(*iterU, nTot);

for(iterV = oldV.begin(); iterV != oldV.end(); iterV++)

{

Vnode addedV = *iterV + Vadd;

//feasibility check

int s1min, s2min, s3min;

minFeasibleValsV(I, J, n, i, s1obs, s2obs, s3obs, s1min, s2min, s3min);

if(addedV.isFeasible(s1obs, s2obs, s3obs, s1min, s2min, s3min))

insertV(&newV, addedV);

}

}//end iterU loop, done defining newV for the i-th primary cluster

//releasing memory for the next i loop

newU.clear();

oldV.clear();

}//end i loop...Vnew is now complete

//finding the number of terms in the conditional distribution

int numt = 0;

for(iterV = newV.begin(); iterV != newV.end(); iterV++)

{

if((*iterV).match(s1obs, s2obs, s3obs))

numt++;

}

//define vecs to hold condtnl dist...note, count can get too big for

unsigned long long

double* t = new double[numt];

double* count = new double[numt];

//filling out the conditional distribution

int ndx=0;

for(iterV = newV.begin(); iterV != newV.end(); iterV++)

{

if((*iterV).match(s1obs, s2obs, s3obs))

{

(*iterV).setDist(t[ndx],count[ndx]); //puts in the value of t,

num for that node

ndx++;

}

}

computePVal(t, count, numt, tobs);

}//end computeDistribution

#include<iostream>

#include<fstream>

#include<cmath>

#include "Unode.h"

using namespace std;

Unode::Unode(int numIn)

{

count=static_cast<double>(numIn);

t = 0;

s1=0;

s2=0;

}

Unode::Unode(double tin, int s1in, int s2in, double numIn)

{

count= numIn;

t=tin;

s1=s1in;

s2=s2in;

}

Unode::Unode(int nijk, int zijk, double xijk)

{

count = choose(nijk, zijk);

t = zijk*xijk;

s1 = zijk;

s2 = zijk*(nijk-zijk);

}

Unode::Unode(const Unode& entry)

{

count = entry.count;

t = entry.t;

s1 = entry.s1;

s2 = entry.s2;

}

bool Unode::operator>(const Unode& entry)

{

if(t>entry.t)

return true;

else if(t==entry.t && s1>entry.s1)

return true;

else if(t==entry.t && s1==entry.s1 && s2>entry.s2)

return true;

return false;

}

bool Unode::operator==(const Unode& entry)

{

return (t==entry.t && s1==entry.s1 && s2==entry.s2);

}

Unode Unode::operator+(const Unode& toAdd)

{

//Note that the correct thing to do with count is multiply, b/c count

//is counting the ways to get a particular outcome.

Unode sum(t + toAdd.t,s1 + toAdd.s1, s2 + toAdd.s2, count*toAdd.count);

return sum;

}

void Unode::print(void)

{

cout << t << " " << s1 << " " << s2 << " " << count << endl;

}

void Unode::augmentCount(const Unode& augmentor)

{

        count += augmentor.count;

}

bool Unode::match(int s1obs, int s2obs)

{

return (s1==s1obs && s2==s2obs);

}

void Unode::setDist(double& tSet, double& countSet)

{

tSet = t;

countSet = count;

}

double Unode::factorial(int number)

{

double temp;

if(number <= 1) return 1;

temp = number * factorial(number - 1);

return temp;

}

double Unode::choose(int top, int bottom)

{

return (factorial(top)/(factorial(bottom)*factorial(top-bottom)));

}

bool Unode::isFeasible(int s1obs, int s2obs)

{

return (s1<=s1obs && s2<=s2obs);

}

void Unode::setVTerms(double& countIn, double& tIn, int& s1In, int& s2In, int&

s3In, int total)

{

countIn = count;

tIn = t;

s1In = s1;

s2In = s2;

s3In = s1*(total-s1); //In TeX shorthand, this is equivalent to \sum_j

z_{ij}(n_{ij}-z_{ij})

}

#include<iostream>

#include<fstream>

#include<cmath>

#include "Vnode.h"

#include "Unode.h"

using namespace std;

Vnode::Vnode(int numIn)

{

count=static_cast<double>(numIn);

t = 0;

s1=0;

s2=0;

s3=0;

}

Vnode::Vnode(double tin, int s1in, int s2in, int s3in, double numIn)

{

count= numIn;

t=tin;

s1=s1in;

s2=s2in;

s3=s3in;

}

Vnode::Vnode(const Vnode& entry)

{

count = entry.count;

t = entry.t;

s1 = entry.s1;

s2 = entry.s2;

s3 = entry.s3;

}

Vnode::Vnode(Unode parent, int total)

{

parent.setVTerms(count, t, s1, s2, s3, total);

}

bool Vnode::operator>(const Vnode& entry)

{

if(t>entry.t)

return true;

else if(t==entry.t && s1>entry.s1)

return true;

else if(t==entry.t && s1==entry.s1 && s2>entry.s2)

return true;

else if(t==entry.t && s1==entry.s1 && s2==entry.s2 && s3>entry.s3)

return true;

return false;

}

bool Vnode::operator==(const Vnode& entry)

{

return (t==entry.t && s1==entry.s1 && s2==entry.s2 && s3==entry.s3);

}

Vnode Vnode::operator+(const Vnode& toAdd)

{

//Note that the correct thing to do with count is multiply, b/c count

//is counting the ways to get a particular outcome.

Vnode sum(t + toAdd.t,s1 + toAdd.s1, s2 + toAdd.s2, s3+ toAdd.s3,

count*toAdd.count);

return sum;

}

void Vnode::print(void)

{

cout << t << " " << s1 << " " << s2 << " " << s3 << " " << count << endl;

}

void Vnode::augmentCount(const Vnode& augmentor)

{

count += augmentor.count;

}

bool Vnode::match(int s1obs, int s2obs, int s3obs)

{

return (s1==s1obs && s2==s2obs && s3==s3obs);

}

void Vnode::setDist(double& tSet, double& countSet)

{

tSet = t;

countSet = count;

}

bool Vnode::isFeasible(int s1obs, int s2obs, int s3obs, int s1min, int s2min,

int s3min)

{

return (s1min<=s1 && s1<=s1obs && s2min<=s2 && s2<=s2obs && s3min<=s3 &&

s3<=s3obs);

}

#ifndef UNODE_H

#define UNODE_H

class Unode

{

//standard variables

double count; //origninally I had this an int, but it can get too big for

unsigned long long's

double t;

int s1;

int s2;

double factorial(int number);

double choose(int top, int bottom);

public:

Unode(int numIn=1); //makes a list with all zeros except num, set to

numIn...use this to make Unew originally

Unode(double tin, int s1in, int s2in, double numIn); //direct constructor

Unode(int nijk, int zijk, double xijk); //constructor from data

Unode(const Unode& entry);

~Unode(){}

//methods for navigating lists

bool operator>(const Unode& entry);

bool operator==(const Unode& entry);

//methods for working with lists in the j-loop

Unode operator+(const Unode& toAdd);

void augmentCount(const Unode& augmentor);

bool match(int s1obs, int s2obs);

void setDist(double& tSet, double& countSet);

bool isFeasible(int s1obs, int s2obs);

//method sets terms; to be used in the V constructor. Total needs to

contain the sum of the n[i][j]'s

//for the ith primary cluster (this needs to also be passed into the

constructor).

void setVTerms(double& countIn, double& tIn, int& s1In, int& s2In, int&

s3In, int total);

void print();

};

#endif

#ifndef VNODE_H

#define VNODE_H

#include "Unode.h"

class Vnode

{

//standard variables

double count; //origninally I had this an int, but it can get too big

for unsigned long long's

double t;

int s1;

int s2;

int s3;

public:

Vnode(int numIn=1); //makes a list with all zeros except num, set to

numIn...use this to make Unew originally

Vnode(Unode parent, int total); //this takes a Unode parent--presumably

from one primary clusters, and makes Vnode

//the argument total needs to contain the sum of n's in that primary cluster.

Vnode(const Vnode& entry);

Vnode(double tin, int s1in, int s2in, int s3in, double numIn); //direct constructor

~Vnode(){}

bool operator>(const Vnode& entry);

bool operator==(const Vnode& entry);

Vnode operator+(const Vnode& toAdd);

void augmentCount(const Vnode& augmentor);

bool match(int s1obs, int s2obs, int s3obs);

void setDist(double& tSet, double& countSet);

bool isFeasible(int s1obs, int s2obs, int s3obs, int s1min, int s2min,

int s3min);

void print();

};

#endif