statease / dexpy

Design of Experiments in Python
https://statease.github.io/dexpy/
Other
27 stars 6 forks source link

2-factor CCD not giving me any star points #64

Closed SandorAlbert closed 1 year ago

SandorAlbert commented 1 year ago

dexpy.ccd.build_ccd() gives me star points when using 3 factors or more, but not when only using 2 factors:

import dexpy.ccd

design = dexpy.ccd.build_ccd(factor_count=2, alpha='orthogonal', center_points=1)
design.index = np.arange(0, len(design))

print(design)

    X1   X2
0 -1.0 -1.0
1 -1.0  1.0
2  1.0 -1.0
3  1.0  1.0
4 -1.0  0.0
5  1.0  0.0
6  0.0 -1.0
7  0.0  1.0
8  0.0  0.0

On the other hand:

import dexpy.ccd

design = dexpy.ccd.build_ccd(factor_count=3, alpha='orthogonal', center_points=1)
design.index = np.arange(0, len(design))

print(design)

          X1        X2        X3
0  -1.000000 -1.000000 -1.000000
1  -1.000000 -1.000000  1.000000
2  -1.000000  1.000000 -1.000000
3  -1.000000  1.000000  1.000000
4   1.000000 -1.000000 -1.000000
5   1.000000 -1.000000  1.000000
6   1.000000  1.000000 -1.000000
7   1.000000  1.000000  1.000000
8  -1.215412  0.000000  0.000000
9   1.215412  0.000000  0.000000
10  0.000000 -1.215412  0.000000
11  0.000000  1.215412  0.000000
12  0.000000  0.000000 -1.215412
13  0.000000  0.000000  1.215412
14  0.000000  0.000000  0.000000

I think there is a miscalculation, since the alpha points should have a | level | of at least 1:

import dexpy.ccd

print(dexpy.ccd.alpha_from_type(factor_count=2, alpha_type="orthogonal", center_points=0))

> 0.9101797211244548

*Edit: I might be totally wrong since this might be the right computation for "orthogonal" designs.

Thanks!

hpanderson commented 1 year ago

The formula for the orthogonal alpha is a function of the number of factorial runs, axial runs, and center points, given by:

$$\alpha = \sqrt[4]{\frac{n_f}{4}\left (\sqrt{n_f+na+n{cp}} - \sqrt{n_f} \right )^2}$$

When you have 1 center point and two factors the correct alpha is indeed 1, and you can see that the quadratic terms are orthogonal to the linear and 2fi terms by examining X'X:

import dexpy.ccd
import numpy as np
import patsy
from dexpy.model import make_quadratic_model

nf = 4
na = 4
nc = 1

alpha = ((nf/4)**.5 * ((nf+na+nc)**.5 - nf**.5)**2)**(.25)
print('alpha', alpha)

design = dexpy.ccd.build_ccd(factor_count=2, alpha='orthogonal', center_points=1)
X = patsy.dmatrix(make_quadratic_model(design.columns), design, return_type="dataframe")
XtX = np.dot(X.T, X)
print(XtX)

[[9. 0. 0. 0. 6. 6.]
 [0. 6. 0. 0. 0. 0.]
 [0. 0. 6. 0. 0. 0.]
 [0. 0. 0. 4. 0. 0.]
 [6. 0. 0. 0. 6. 4.]
 [6. 0. 0. 0. 4. 6.]]

Incidentally, this is also the D-optimal design for 2 factors in 9 runs.

If you perform the same exercise with 0 center points the alpha does go down to ~0.91 which is interesting, but appears to be correct.