def analyze_theory(V, x0list=[], plot=False):
    """ Extract ground-state energy E0 and psi**2 for potential V. """
    # initialize path integral
    T = 4.
    ndT = 8.         # use larger ndT to reduce discretization error (goes like 1/ndT**2)
    neval = 3e5   # should probably use more evaluations (10x?)
    nitn = 6
    alpha = 0.1     # damp adaptation

    # create integrator and train it (no x0list)
    integrand = PathIntegrand(V=V, T=T, ndT=ndT)
    integ = vegas.Integrator(integrand.region, alpha=alpha)
    integ(integrand, neval=neval, nitn=nitn / 2, alpha=2 * alpha)

    # evaluate path integral with trained integrator and x0list
    integrand = PathIntegrand(V=V, x0list=x0list, T=T, ndT=ndT)
    results = integ(integrand, neval=neval, nitn=nitn, alpha=alpha)
    print(results.summary())
    E0 = -np.log(results['exp(-E0*T)']) / T
    print('Ground-state energy = %s    Q = %.2f\n' % (E0, results.Q))

    if len(x0list) <= 0:
        return E0
    psi2 = results['exp(-E0*T) * psi(x0)**2'] / results['exp(-E0*T)']
    print('%5s  %-12s %-10s' % ('x', 'psi**2', 'sho-exact'))
    print(27 * '-')
    for i, (x0i, psi2i) in enumerate(zip(x0list, psi2)):
        exact = np.exp(- x0i ** 2) / np.sqrt(np.pi) #* np.exp(-T / 2.)
        print(
            "%5.1f  %-12s %-10.5f"
            % (x0i, psi2i, exact)
            )
    if plot:
        plot_results(E0, x0list, psi2, T)
    return E0
Extract ground-state energy E0 and psi**2 for potential V.