def orth_chol(order, dist, normed=True, sort="GR", cross_truncation=1., **kws):
    dim = len(dist)
    basis = chaospy.poly.basis(
        start=1, stop=order, dim=dim, sort=sort,
        cross_truncation=cross_truncation,
    )
    length = len(basis)
    cholmat = chaospy.chol.gill_king(chaospy.descriptives.Cov(basis, dist))
    cholmat_inv = numpy.linalg.inv(cholmat.T).T
    if not normed:
        diag_mesh = numpy.repeat(numpy.diag(cholmat_inv), len(cholmat_inv))
        cholmat_inv /= diag_mesh.reshape(cholmat_inv.shape)
    coefs = numpy.empty((length+1, length+1))
    coefs[1:, 1:] = cholmat_inv
    coefs[0, 0] = 1
    coefs[0, 1:] = 0
    expected = -numpy.sum(
        cholmat_inv*chaospy.descriptives.E(basis, dist, **kws), -1)
    coefs[1:, 0] = expected
    coefs = coefs.T
    out = {}
    out[(0,)*dim] = coefs[0]
    for idx in range(length):
        index = basis[idx].keys[0]
        out[index] = coefs[idx+1]
    polynomials = chaospy.poly.Poly(out, dim, coefs.shape[1:], float)
    return polynomials