diff --git a/PLAN.md b/PLAN.md
index 4464e63..3a77b99 100644
--- a/PLAN.md
+++ b/PLAN.md
@@ -94,12 +94,19 @@ replace project documentation.
its original authorship, date, Keister example, code, mathematics, output,
references, and figure.
- [x] Establish listing metadata and a discoverable Blog index.
-- [ ] Inventory the existing 18 posts and their assets without modifying the
+- [x] Inventory the existing 18 posts and their assets without modifying the
source repository.
-- [ ] Define Quarto post metadata, URL preservation, authorship, categories,
- and redirect requirements.
-- [ ] Convert and review posts in small batches.
-- [ ] Verify dates, authors, links, code, mathematics, images, and legacy URLs.
+- [x] Define consistent Quarto post metadata, paths, authorship, and categories.
+- [x] Use one `date` value per post as the last-revised date, label it “Last
+ revised,” and sort the archive by that value without separate first-published
+ metadata.
+- [x] Convert the complete 18-post archive, including MkDocs snippets, callouts,
+ image groups, internal links, and mathematical delimiters.
+- [x] Preserve the original visible content of 17 posts and retain the approved
+ living-article revisions to “Why Add Q to MC?”.
+- [x] Verify last-revised dates, authors, local links, code, mathematics, and all
+ 48 images in a complete local render and browser review.
+- [ ] Define and verify redirects or canonical handling for legacy URLs.
- [ ] Publish the complete archive only after collaborator review.
### Phase 6 — Build news and events
@@ -167,4 +174,5 @@ replace project documentation.
with QMCSoftware collaborators.
- [ ] Record requested changes and resolve launch-blocking issues.
- [ ] Agree on ownership and cadence for blog, news, and community updates.
-- [ ] Approve the prototype before beginning the 18-post blog migration.
+- [x] Authorize and locally validate the 18-post blog migration.
+- [ ] Approve the completed local archive before publishing it for review.
diff --git a/STATUS.md b/STATUS.md
index c91158b..048dbcd 100644
--- a/STATUS.md
+++ b/STATUS.md
@@ -32,19 +32,30 @@
- [x] Added repository guidance, roadmap, ignore rules, and local instructions.
- [x] Rendered the complete 13-page site and verified local link targets and
rendered `CNAME` preservation.
+- [x] Migrated the complete 18-post QMCPy blog archive into self-contained
+ Quarto posts while retaining the technical documentation and notebooks in
+ the QMCPy repository.
+- [x] Standardized blog chronology on one last-revised `date` per post, with no
+ separate first-publication field, and retained the approved revised version
+ of “Why Add Q to MC?”.
+- [x] Preserved and published all 48 blog images in the local render, expanded
+ MkDocs code snippets, converted callouts and image groups, and normalized
+ legacy mathematical delimiters for Quarto.
+- [x] Rendered all 30 site pages and completed local link, desktop, mobile, and
+ per-post browser checks for the 18-post archive.
## Underway
- [ ] Collaborator review of the information architecture, visual direction,
and provisional homepage copy.
-- [ ] Browser-based responsive and accessibility review.
+- [ ] Accessibility review beyond the completed desktop and mobile browser
+ checks.
- [ ] Live GitHub Pages workflow and custom-domain verification.
## Next recommended task
-- [ ] Review the rendered umbrella prototype with QMCSoftware collaborators and
- approve the visual direction and public copy before adding content or
- migrating posts.
+- [ ] Review and approve the completed local 18-post archive before any branch
+ is pushed or a pull request is opened.
## Decisions and scope boundaries
@@ -53,8 +64,13 @@
link to their respective repositories and documentation.
- The site uses `qmcsoftware.org` as its canonical URL and custom domain;
redirect behavior from the former domain still requires verification.
-- Two earlier blog posts are migrated as a pattern; the remainder of the
- archive is outside the current prototype.
+- The complete 18-post QMCPy blog archive is migrated locally; publishing it
+ remains approval-gated.
+- Blog cards and title blocks use each post's last-revised date; the archive is
+ sorted newest revision first and does not carry a separate first-published
+ date.
+- “Why Add Q to MC?” remains the approved living article; the other 17 posts
+ preserve the visible content of their source articles.
- No Reveal.js project, speculative infrastructure, or Python dependency file
is included.
- Generated `_site/` output remains outside the `main` branch.
diff --git a/_quarto.yml b/_quarto.yml
index 2f78e4b..b3fa00b 100644
--- a/_quarto.yml
+++ b/_quarto.yml
@@ -12,6 +12,7 @@ project:
- CNAME
- assets/**
- data/**
+ - blog/**/figures/**
website:
title: "QMCSoftware"
diff --git a/blog/_metadata.yml b/blog/_metadata.yml
new file mode 100644
index 0000000..bbebb9b
--- /dev/null
+++ b/blog/_metadata.yml
@@ -0,0 +1,2 @@
+language:
+ title-block-published: "Last revised"
diff --git a/blog/accelerating-qmcpy-notebook-tests-with-parsl/figures/parsl-four-worker-speedup.png b/blog/accelerating-qmcpy-notebook-tests-with-parsl/figures/parsl-four-worker-speedup.png
new file mode 100644
index 0000000..3c6cd4e
Binary files /dev/null and b/blog/accelerating-qmcpy-notebook-tests-with-parsl/figures/parsl-four-worker-speedup.png differ
diff --git a/blog/accelerating-qmcpy-notebook-tests-with-parsl/figures/parsl-worker-scaling.png b/blog/accelerating-qmcpy-notebook-tests-with-parsl/figures/parsl-worker-scaling.png
new file mode 100644
index 0000000..63e3b1e
Binary files /dev/null and b/blog/accelerating-qmcpy-notebook-tests-with-parsl/figures/parsl-worker-scaling.png differ
diff --git a/blog/accelerating-qmcpy-notebook-tests-with-parsl/index.qmd b/blog/accelerating-qmcpy-notebook-tests-with-parsl/index.qmd
new file mode 100644
index 0000000..8f4e97a
--- /dev/null
+++ b/blog/accelerating-qmcpy-notebook-tests-with-parsl/index.qmd
@@ -0,0 +1,101 @@
+---
+title: "Accelerating QMCPy Notebook Tests with Parsl"
+author: "Joshua Jay Herman, Brandon Sharp, and Sou-Cheng Choi"
+date: 2025-11-29
+date-format: "MMMM D, YYYY"
+description: "How Parsl parallelism accelerates QMCPy notebook testing while preserving the existing Testbook workflow."
+categories:
+ - "QMCPy"
+ - "Testing"
+ - "Parallel Computing"
+image: figures/parsl-worker-scaling.png
+---
+
+This post describes how Parsl parallelism can accelerate QMCPy notebook testing while preserving the structure of the existing testbook workflow.
+
+## Introduction
+
+Notebook regression testing ensures that interactive examples and analyses remain correct and reproducible, catching regressions introduced by changes in code, dependencies, or execution environments. For QMCPy [1], this process is both massively parallel and resource-intensive because of the number and complexity of its notebooks.
+
+This blog post summarizes our work on accelerating notebook regression testing using Testbook-based tests, which can be viewed as notebook-level unit tests. The work was presented in our ParslFest 2025 talk [2], and this post outlines directions for further development.
+
+The presentation slides are available as [Parsl Testbook Speedup](https://qmcsoftware.github.io/QMCSoftware/demos/talk_paper_demos/Parslfest_2025/Parsl_Testbook_Speedup.pdf).
+
+## Methodology
+
+Our choice to adopt Testbook [3] is motivated by its ability to execute Jupyter notebooks directly within a test environment, enabling fine-grained validation of both code cells and notebook state. Testbook also integrates cleanly with our existing testing directory structure, where other unit tests are organized without requiring full notebook execution. This preserves modularity, simplifies debugging, and avoids unnecessary duplication of logic.
+
+To support scalable notebook testing, we developed a lightweight yet flexible test harness that enables Parsl [4] to orchestrate Testbook-based unit tests. By treating each notebook test as an independent Parsl app, the harness realizes an embarrassingly parallel workflow suitable for local multiprocessing, HPC schedulers, or cloud environments.
+
+The harness coordinates three primary components to achieve reproducible, high-throughput notebook testing:
+
+- **Continuous Integration (CI):** A GitHub Actions workflow prepares the
+execution environment, including Conda environment creation, minimal LaTeX installation, optional swap configuration, project dependency installation, and test-target execution. This ensures consistent, version-controlled execution across platforms.
+- **Parsl controller and workers:** Parsl provisions local or remote executors,
+including processes, threads, or cluster jobs, and schedules notebook tests as independent tasks. This enables parallel execution with configurable concurrency limits, resource profiles, and executor backends.
+- **Testbook runner and artifact collection:** Each worker executes its
+assigned notebook tests through Testbook. Outputs, execution logs, error traces, generated figures, and notebook artifacts with executed cells are returned to the Parsl controller and can be uploaded by CI for inspection, provenance tracking, and debugging.
+
+Key features of the harness include pinned Conda environment specifications for reproducibility, customizable Parsl executors, timeout and retry policies for handling flaky or long-running tests, and centralized logging to streamline diagnosis of failures. Together, these components provide a robust framework for scalable, automated validation of computational notebooks.
+
+## Results
+
+To establish a performance baseline, we first measured the wall-clock time required to execute a representative subset of demo notebooks sequentially. After extending test coverage to include syntax-validation checks and additional notebooks, we repeated the experiment under the parallel Testbook-Parsl workflow. Across these configurations, we observed a consistent 3.0-fold speedup, demonstrating that notebook-based tests parallelize cleanly and benefit substantially from concurrent execution.
+
+
+
+ Figure 1: Speedup achieved by running Testbook-based notebook tests under Parsl with four workers compared to sequential execution.
+
+
+We also measured how execution time and speedup changed as the number of Parsl workers increased. The gains improved substantially from two to four workers, then flattened as overheads and resource limits began to dominate.
+
+
+
+ Figure 2: Speedup achieved by running Testbook-based notebook tests under Parsl with different worker counts compared to sequential execution.
+
+
+All tests were executed on a Linux system with AMD64 architecture and 16 CPU cores. When run in continuous integration, the workflow executes the same test suite on a GitHub Actions `ubuntu-latest` runner. Users may reproduce the parallel Testbook workflow locally by running:
+
+```bash
+make booktests_parallel_no_docker
+```
+
+The experiment covered the following notebook-test scope:
+
+| Item | Count | Notes |
+|---|---:|---|
+| Demo notebooks in the experiment | 33 | Representative notebooks from `demos/`, including subfolders |
+| Generated test files in `test/booktests/tb_*.py` | 32 | Each `tb_*.py` tests a single demo notebook |
+| Notebooks executed by the booktests workflow | 32 | Workflow runs the `make` targets and executes one generated test file per notebook |
+
+## Runner Configuration Notes
+
+The GitHub `ubuntu-latest` runner typically provides 2 virtual CPUs per job. The workflow allocates a 12 GB swap file to mitigate transient memory spikes during notebook execution and reduce the likelihood of out-of-memory failures.
+
+## CI Tests
+
+The continuous integration workflow automates the execution of notebook-based tests and prepares a controlled environment for reproducible runs. The workflow checks out the repository, sets up Miniconda, installs the project's test extras with `pip install -e .[test]` and optional components such as `test_torch`, `test_gpytorch`, `test_botorch`, and `test_umbridge`, and creates a 12 GB swap file early in the job to reduce out-of-memory failures for notebooks with heavy memory demands.
+
+To ensure notebook tests are present and up to date, the workflow invokes `make check_booktests` and `make generate_booktests`, which confirm or regenerate `test/booktests/tb_*.py` files. The final step triggers either `make booktests_parallel_no_docker` for parallel execution or `make booktests_no_docker` for sequential execution. These targets run the generated tests via `pytest` and Testbook.
+
+For diagnostics, setting `ACTIONS_STEP_DEBUG=true` increases log verbosity. Timing and memory usage can be captured by wrapping the `make` call with `/usr/bin/time -v` and uploading the resulting logs using the `actions/upload-artifact@v4` step.
+
+Parallel notebook execution must respect the resource limits of GitHub runners, which typically provide only two CPUs. When notebooks request more workers, adapting the code to use `max_workers = min(8, os.cpu_count() or 1)` helps ensure compatibility. In cases of memory pressure, the workflow may fall back to the sequential `booktests_no_docker` target or reduce parallelism inside the `Makefile`.
+
+The workflow is readily extensible: caching package installations with `actions/cache` accelerates subsequent runs; notebook outputs and HTML artifacts can be uploaded for failure analysis; and a simple CSV timing log in `test/booktests/` can be collected to track notebook performance over time.
+
+## Further Work
+
+These results indicate that we should extend this style of parallel testing to doctests and `pytest` tests. Because many developers have multicore processors, parallel local testing can improve individual productivity while still demonstrating that no regressions have been introduced.
+
+Feedback from ParslFest participants also highlighted that the system is quite general. This suggests that a distributed test system could benefit Parsl users by enabling them to distribute their own test workloads. Future work could expand this approach to Python doctests and unit testing with `pytest` or `unittest`, in addition to testing Jupyter notebooks.
+
+## References
+
+1. Choi, S.-C. T., Hickernell, F. J., McCourt, M., Rathinavel, J., & Sorokin,
+ A. QMCPy: A quasi-Monte Carlo Python Library. [https://qmcpy.org](https://qmcpy.org).
+2. ParslFest 2025 presentation materials:
+ [Parsl Testbook Speedup](https://qmcsoftware.github.io/QMCSoftware/demos/talk_paper_demos/Parslfest_2025/Parsl_Testbook_Speedup.pdf).
+3. Testbook documentation. [https://testbook.readthedocs.io/](https://testbook.readthedocs.io/).
+4. Babuji, Y. et al. Parsl: Pervasive Parallel Programming in Python.
+*Proceedings of the 28th International Symposium on High-Performance Parallel and Distributed Computing* (2019).
diff --git a/blog/accelerating-rare-event-reliability-simulations-for-cerns-large-hadron-collider-using-qmcpy/figures/results.png b/blog/accelerating-rare-event-reliability-simulations-for-cerns-large-hadron-collider-using-qmcpy/figures/results.png
new file mode 100644
index 0000000..5a7704b
Binary files /dev/null and b/blog/accelerating-rare-event-reliability-simulations-for-cerns-large-hadron-collider-using-qmcpy/figures/results.png differ
diff --git a/blog/accelerating-rare-event-reliability-simulations-for-cerns-large-hadron-collider-using-qmcpy/figures/system-tree.png b/blog/accelerating-rare-event-reliability-simulations-for-cerns-large-hadron-collider-using-qmcpy/figures/system-tree.png
new file mode 100644
index 0000000..e83a582
Binary files /dev/null and b/blog/accelerating-rare-event-reliability-simulations-for-cerns-large-hadron-collider-using-qmcpy/figures/system-tree.png differ
diff --git a/blog/accelerating-rare-event-reliability-simulations-for-cerns-large-hadron-collider-using-qmcpy/figures/timeline.png b/blog/accelerating-rare-event-reliability-simulations-for-cerns-large-hadron-collider-using-qmcpy/figures/timeline.png
new file mode 100644
index 0000000..0819620
Binary files /dev/null and b/blog/accelerating-rare-event-reliability-simulations-for-cerns-large-hadron-collider-using-qmcpy/figures/timeline.png differ
diff --git a/blog/accelerating-rare-event-reliability-simulations-for-cerns-large-hadron-collider-using-qmcpy/index.qmd b/blog/accelerating-rare-event-reliability-simulations-for-cerns-large-hadron-collider-using-qmcpy/index.qmd
new file mode 100644
index 0000000..4b91e0d
--- /dev/null
+++ b/blog/accelerating-rare-event-reliability-simulations-for-cerns-large-hadron-collider-using-qmcpy/index.qmd
@@ -0,0 +1,80 @@
+---
+title: "Accelerating Rare-event Reliability Simulations for CERN's Large Hadron Collider using QMCPy"
+author: "Milosz Blaszkiewicz"
+date: 2023-08-26
+date-format: "MMMM D, YYYY"
+description: "Using QMCPy to accelerate rare-event Monte Carlo reliability simulations in CERN's AvailSim4 framework."
+categories:
+ - "QMCPy"
+ - "Reliability"
+ - "Quasi-Monte Carlo"
+image: figures/results.png
+---
+
+This post shows how QMCPy can improve rare-event Monte Carlo reliability simulations in CERN's AvailSim4 framework.
+
+In this blog post, we share an example for using the QMCPy package to accelerate rare-event Monte Carlo (MC) simulations in AvailSim4 [1]. The effort is part of a more general study of advanced MC methods for reliability studies of the CERN Machine Protection group.
+
+## Introduction
+
+The European Organization for Nuclear Research, CERN, is home to the largest particle accelerator in the world, the Large Hadron Collider (LHC). The machine is producing and recording data from high-energy collisions of proton beams, allowing scientists all over the world to test theoretical models and hypotheses. The unprecedented beam energies create potential for pushing further the very boundaries of human knowledge and answering the most fundamental questions regarding the origins of the universe.
+
+The LHC consists of many sophisticated systems with responsibilities such as injecting the particles into the beam orbit, maintaining beams on very precise tracks, or cooling down the superconducting magnets and radio frequency cavities, to name just a few. The Machine Protection group is primarily focused on two types of systems: those protecting the superconducting magnets and their circuits and those protecting the accelerator equipment from damage due to the circulating high-energy beams. These systems are critical as their failures may cause severe damage to the accelerator.
+
+Reliability engineering offers a wide range of methods to quantify potential risks and their consequences. Probabilistic methods such as Monte Carlo (MC) simulations are used to assess risks of complex systems for which no analytic solution can be derived. However, MC simulation may come at a significant computational cost. In this short post, we present the AvailSim4 tool in which we have implemented a Quasi-Monte Carlo extension to make the computations more efficient. The open-source implementation of the framework utilizes the QMCPy package.
+
+## AvailSim4
+
+AvailSim4 is an open-source framework developed in the Machine Protection group of CERN's Technology Department. AvailSim4 provides an environment for availability simulations with several features specifically developed for particle accelerator applications.
+
+
+
+ Simplified AvailSim4 system model with compound and basic components.
+
+
+An AvailSim4 model can include several sub-systems, each described by a failure probability distribution, recovery distribution, and functional dependencies. The structure of those dependencies forms a tree, which models the entire system. A simplified example of such a system is presented in the picture above. The overall component **System** consists of two children: basic **A** without children and compound **B**, which further splits into two components, **B1** and **B2**. The relation between the two is defined to be guided by **OR** logic: component **B** is working if either of its two children is operational.
+
+
+
+ Visualization of a simulated event timeline and supporting component timelines.
+
+
+AvailSim4 combines MC with a Discrete Event Simulation (DES) approach. This means that each MC iteration is a realization of a random timeline of events. A visualization of a timeline, with supporting timelines of individual components, is displayed in the plot above. The overall analysis is the result of how often and for how long individual components fail. In this sense, the framework is a standard example of an MC simulation: it runs multiple instances of the experiment with different random numbers and obtains estimations of quantities of interest by calculating their averages.
+
+The main difficulty is the shift to a rare-event regime, where events of interest occur in only a small fraction of iterations. Computation efficiency can be improved on two separate levels: code optimizations, such as profiling, distributed computing, and precompiling critical functions, and algorithm enhancements. Applying Quasi-Monte Carlo is a member of the latter group.
+
+## Challenges of Quasi-Monte Carlo in AvailSim4
+
+Using Quasi-Monte Carlo did not require substantial changes in the Monte Carlo implementation in the case of AvailSim4. The bigger change is on the conceptual level. Understanding the differences is essential for interpretation of our results.
+
+However, the use of QMC methods comes with a dimensionality limitation, which is important in the studied use case. First, the computational efficiency improvement comes from using samples covering the problem space more evenly. In our use case, that space should be viewed in the MC sense rather than a single DES sense. This means that consecutive samples will be employed across different DES iterations, and the need for more random values in an individual iteration will be addressed by generating samples of multiple dimensions. This aspect will be further discussed below.
+
+## Results
+
+The test case is made of a very simple system, consisting of a few redundant components that have the same properties and fail at times drawn from an exponential probability distribution. The aspect that changes between the three presented cases is the number of those components. In such a scenario, the more components are needed, the less likely a critical failure is.
+
+
+
+ Accuracy progression and execution time comparison for MC and QMC modes across rare-event test cases.
+
+
+In the left-hand side plots above, we see the progression of the accuracy as the number of DES iterations increases. Orange lines represent QMC results, while blue ones are results of the MC mode. In only one case, 5 components and 100 iterations, the QMC mode is less accurate. All remaining test cases show that the orange line is closer to the value to which both lines eventually converge. The execution time comparison is featured on the right side. Results are also relatively stable: QMC adds a small overhead at the beginning to generate a large matrix of random numbers before iteration rather than in it, however it is visible only in the case of small numbers of iterations. This overhead does not eliminate the advantage of the method, which is coming from using fewer iterations required to obtain certain accuracy. Also, the more iterations are completed, the smaller the relative difference, as generating random values takes place only once.
+
+All results presented in this section are further discussed in [2]. This includes additional test cases and a comparison of QMC with Importance Splitting, another method to significantly speed up MC simulations.
+
+Another aspect is limitations of the approach and using QMC in general. It has already been said that samples are multidimensional, so that each one contains enough random values for all components in the DES. However, there is an additional complication. A significant element of all availability and reliability simulations is that components may be repaired and returned to their fully operational state an indefinite number of times. Assigning each component only a single failure time is a solution that falls short in those terms. Instead, the existing implementation provisions more random numbers, by assigning more dimensions, for each component. Whenever a given element fails, the next failure time is taken as a value of the next dimension assigned to it.
+
+This fact brings the most serious limitation of the approach. The number of available random failure time values needs to be decided prior to commencing simulations and must assume the worst case, so that no component runs out of failure times before finishing its lifetime. When some of the components fail relatively often, the total number of dimensions will often end up close to the current limits of low-discrepancy sequence generators. This also adds to the overhead at the beginning: we need to pre-emptively generate many more values than when generating them only when needed, i.e., as it is done with a standard pseudo-random number generator.
+
+## Conclusions
+
+During this study, the QMC method evolved from a proof-of-concept addition to a fully implemented feature of AvailSim4. The most important advantage of QMC methods is that the change from crude MC is almost transparent for the users: no further information or inputs are required. It would be fully transparent had there not been the limitation caused by the number of dimensions of each sample, which is something to which users need to pay attention.
+
+In our tests of the rare-events scenarios, the gains are visible. The results' accuracy increased as their variance diminished. However, getting orders of magnitude improvements is strictly impossible, as the method of obtaining results is still based on simple calculation of averages, such as the numbers of events of interest.
+
+## References
+
+1. AvailSim4 GitLab repository.
+ [https://gitlab.cern.ch/availsim4/availsim4](https://gitlab.cern.ch/availsim4/availsim4).
+2. Blaszkiewicz, M. *Methods to optimize rare-event Monte Carlo
+reliability simulations for Large Hadron Collider Protection Systems*. MSc Thesis. [https://cds.cern.ch/record/2808520](https://cds.cern.ch/record/2808520).
diff --git a/blog/analysis-of-qmc-efficiency-for-asian-option-pricing/figures/qmc_points_required_combined.png b/blog/analysis-of-qmc-efficiency-for-asian-option-pricing/figures/qmc_points_required_combined.png
new file mode 100644
index 0000000..c32731e
Binary files /dev/null and b/blog/analysis-of-qmc-efficiency-for-asian-option-pricing/figures/qmc_points_required_combined.png differ
diff --git a/blog/analysis-of-qmc-efficiency-for-asian-option-pricing/figures/qmc_time_required_combined.png b/blog/analysis-of-qmc-efficiency-for-asian-option-pricing/figures/qmc_time_required_combined.png
new file mode 100644
index 0000000..df5e9ce
Binary files /dev/null and b/blog/analysis-of-qmc-efficiency-for-asian-option-pricing/figures/qmc_time_required_combined.png differ
diff --git a/blog/analysis-of-qmc-efficiency-for-asian-option-pricing/index.qmd b/blog/analysis-of-qmc-efficiency-for-asian-option-pricing/index.qmd
new file mode 100644
index 0000000..5094bc3
--- /dev/null
+++ b/blog/analysis-of-qmc-efficiency-for-asian-option-pricing/index.qmd
@@ -0,0 +1,151 @@
+---
+title: "Analysis of Quasi-Monte Carlo Efficiency for Asian Option Pricing"
+author: "Karm Dave"
+date: 2025-07-16
+date-format: "MMMM D, YYYY"
+description: "A comparison of IID Monte Carlo and QMCPy QMC methods for arithmetic Asian option pricing."
+categories:
+ - "QMCPy"
+ - "Finance"
+ - "Quasi-Monte Carlo"
+image: figures/qmc_points_required_combined.png
+---
+
+This post compares IID Monte Carlo and QMCPy QMC methods for arithmetic Asian option pricing, focusing on sample counts, timings, and empirical stability.
+
+## Introduction: The Asian Option Pricing Challenge
+
+Asian options are especially difficult to price because their payoffs involve the average price of an underlying asset over a period of time. In comparison to simpler European options, arithmetic (or average) Asian options generally do not have a simple closed-form solution because of their path dependence. As a result, they are commonly approximated using numerical methods, including Monte Carlo (MC) simulation.
+
+In standard MC simulations, the model simulates a large number of asset price paths, where each path is generated by a different set of pseudorandom numbers, and averages the resulting discounted payoffs. The root-mean-square error of standard MC is typically $O(N^{-1/2})$, as derived from the Central Limit Theorem. This convergence rate is dimension-independent, but it can require many samples when high accuracy is needed.
+
+Quasi-Monte Carlo (QMC) methods provide an alternative to MC methods by substituting deterministic low-discrepancy sequences (e.g., Sobol’, Halton) or point sets (e.g., lattice rules) in place of pseudorandom numbers. These constructions cover the unit cube more evenly than independent pseudorandom points. In favorable settings, especially when the integrand is sufficiently smooth and has moderate effective dimension, this can produce faster convergence than standard MC. In this work, we compare MC and QMC methods for pricing an arithmetic Asian call option and report median performance across repeated runs.
+
+## Mathematical Framework
+
+We consider an arithmetic Asian call option with strike $K$ and expiry
+$T$. The underlying asset price $S(t)$ is observed at $d$ discrete
+times $t_j = jT/d$. We use the trapezoidal arithmetic mean convention,
+so the discounted payoff is:
+
+$$
+\text{Payoff} =
+\max\left(\frac{S(0)/2 + \sum_{j=1}^{d-1} S(t_j) + S(T)/2}{d} - K, 0\right)e^{-rT}
+$$
+
+where $r$ is the risk-free rate. The goal is to estimate the expected value of this payoff.
+
+Asset prices are modeled using Geometric Brownian Motion (GBM) under the risk-neutral measure:
+
+$$
+dS(t) = rS(t)dt + \sigma S(t)dW(t)
+$$
+
+where $\sigma$ is volatility and $W(t)$ is a Wiener process.
+
+## Simulation Setup
+
+We used [QMCPy](https://qmcpy.org) to compare standard MC (IID) against QMC (Sobol’, Halton, Lattice) for pricing an Asian call option with parameters:
+
+- $S_0 = \$120$
+- $K = \$130$
+- $T = 1$ year
+- $r = 0.02$
+- $\sigma = 0.50$
+- $d = 12$ monthly observations
+
+The integration dimension is $d = 12$.
+
+We used adaptive algorithms (`CubMCG`, `CubQMCSobolG`, `CubQMCCLT`, `CubQMCLatticeG`) to achieve target absolute error tolerances $\epsilon$ from \$0.50 to as low as \$0.0005. For each method and tolerance, we ran the simulation 25 times to collect data on sample needs and processing time. The numbers we report for sample counts and timings show the median value across these runs. We include the 25th and 75th percentiles in graphs to show the variability in results.
+
+## QMCPy Sobol’ Example
+
+The following example shows how to use Sobol’ sequence in QMCPy for an Asian option pricing simulation.
+
+```python
+import qmcpy as qp
+
+d = 12
+
+# Point Generators (re-instantiated for each run)
+SobolPoints = qp.Sobol(dimension=d, seed=7)
+IIDPoints = qp.IIDStdUniform(dimension=d) # etc.
+
+# Define Asian Option Measure (Sobol’)
+ArithMeanCallSobol = qp.AsianOption(
+ SobolPoints,
+ volatility=0.50,
+ start_price=120,
+ strike_price=130,
+ interest_rate=0.02,
+ t_final=1,
+ call_put='call',
+ asian_mean="ARITHMETIC",
+)
+
+# Define similarly for IID, Halton, Lattice...
+
+# Run Integration (example for one run, one tolerance)
+price, data = qp.CubQMCSobolG(
+ ArithMeanCallSobol,
+ abs_tol=0.01,
+ rel_tol=0,
+).integrate()
+
+print(f"Sobol’ Price: {price:.4f}, Samples: {data.n_total}")
+
+# Example output: Sobol’ Price: 10.3167, Samples: 16384
+# Full script iterates this over tolerances and multiple runs.
+```
+
+## Results: QMC Efficiency Demonstrated
+
+As Table **1** demonstrates, different method and stopping-criterion combinations required different numbers of samples to achieve the desired accuracy. QMC stopping criteria often add samples in powers of two, so Sobol’ and lattice runs can have stable sample counts for a given tolerance. These medians show the central tendency across a range of results for IID (MC) and Halton (with `CubQMCCLT`). The comparison should therefore be read as a comparison of the full QMCPy integration workflows, not only as a comparison of point generators. Sobol’ and Lattice perform better than Halton in part because `CubQMCSobolG` and `CubQMCLatticeG` use decay estimates of Walsh or Fourier coefficients to select sample sizes, whereas `CubQMCCLT` for Halton uses multiple randomizations and a central-limit-theorem-based stopping criterion.
+
+The following table summarizes the median samples required for each target absolute error tolerance. `N/A` means that the tolerance was not met practically by IID in this experiment.
+
+| Error $\epsilon$ | Sobol’ | IID (MC) | Halton | Lattice |
+|---:|---:|---:|---:|---:|
+| 0.5000 | 1,024 | 38,602 | 4,096 | 1,024 |
+| 0.3053 | 1,024 | 85,400 | 4,096 | 1,024 |
+| 0.1864 | 1,024 | 185,561 | 4,096 | 1,024 |
+| 0.1138 | 1,024 | 474,908 | 4,096 | 1,024 |
+| 0.0695 | 2,048 | 1,265,302 | 8,192 | 2,048 |
+| 0.0424 | 4,096 | N/A | 16,384 | 4,096 |
+| 0.0259 | 8,192 | N/A | 32,768 | 8,192 |
+| 0.0158 | 8,192 | N/A | 32,768 | 8,192 |
+| 0.0097 | 16,384 | N/A | 65,536 | 16,384 |
+| 0.0059 | 32,768 | N/A | 131,072 | 32,768 |
+| 0.0036 | 65,536 | N/A | 262,144 | 65,536 |
+| 0.0022 | 65,536 | N/A | 262,144 | 65,536 |
+| 0.0013 | 131,072 | N/A | 524,288 | 131,072 |
+| 0.0008 | 262,144 | N/A | 1,048,576 | 262,144 |
+| 0.0005 | 524,288 | N/A | 1,048,576 | 524,288 |
+
+The final high-precision estimates at $\epsilon = 0.0005$ from the last simulation run were consistent across QMC methods: Sobol’ (\$10.3173), Lattice (\$10.3173), Halton (\$10.3172). Standard MC could not produce reliable estimates below about $\epsilon = 0.05$ in the last simulation run, indicating practical limitations at higher precision levels.
+
+Figure 1 visually underscores the efficiency gap, plotting median requirements with 25th-75th percentile ranges to indicate variability.
+
+
+
+
+
+
+
+
+Median sample and time requirements versus error tolerance, with 25th-75th percentile ranges shown as shaded areas based on 25 runs. Left: median samples required. Right: median runtime.
+
+
+Several observations can be drawn from the results. For the tolerance $\epsilon=0.0695$, both Sobol' and Lattice required substantially fewer samples than IID MC, using roughly 618 times fewer paths in the median comparison ($N_{\text{Sobol'}} = 2,048$ vs. $N_{\text{IID}} \approx 1,265,302$). At tighter tolerances, the gap in sample requirements becomes larger in this experiment. Sobol' and Lattice produced identical median sample counts across the tested tolerances, which suggests similar practical behavior for this particular option pricing setup and these stopping criteria. Halton performed better than MC, but required more samples than Sobol' and Lattice at the tested tolerances; this partly reflects the use of the `CubQMCCLT` stopping criterion with multiple randomizations. Standard MC did not reach the tighter tolerances in the last simulation run under the tested computational budget. Overall, the empirical results are consistent with the expected advantage of low-discrepancy constructions for sufficiently regular integrands. The variability plots also show that the QMC workflows had more stable sample counts than MC in this experiment.
+
+## Performance Edge and Practical Implications
+
+The better sample efficiency of QMC, evidenced by median performance across 25 runs, has practical implications for pricing and risk analysis. For this example, reducing the number of simulated paths also reduced runtime, making tighter tolerances more feasible than with the tested MC workflow. This matters when the payoff evaluation is expensive, when many parameter settings must be evaluated, or when computations are run repeatedly in production or cloud environments. The cost of generating QMC points is small relative to the cost of simulating paths and evaluating payoffs in this experiment, so the lower sample counts translate directly into computational savings.
+
+## Conclusion
+
+For this arithmetic Asian option example, the tested QMC workflows achieved the target tolerances with far fewer samples than the tested standard MC workflow. Sobol' sequences and Lattice rules were especially effective in this setup, while Halton also improved on MC but required more samples under the `CubQMCCLT` stopping criterion. These results do not imply that QMC will dominate MC for every path-dependent derivative, but they show why low-discrepancy sampling is often attractive when the integrand is sufficiently regular and the effective dimension is manageable. Repeated runs using medians and quantiles provide a more stable view of the comparison than a single run alone.
+
+## References
+
+- Choi, S.-C. T., et al., “QMCPy: A Quasi-Monte Carlo Python Library,” [https://qmcpy.org](https://qmcpy.org).
diff --git a/blog/bayesian-stopping-criteria/index.qmd b/blog/bayesian-stopping-criteria/index.qmd
new file mode 100644
index 0000000..b090bf0
--- /dev/null
+++ b/blog/bayesian-stopping-criteria/index.qmd
@@ -0,0 +1,261 @@
+---
+title: "Bayesian Stopping Criteria"
+author: "Jagadeeswaran Rathinavel and Fred J. Hickernell"
+date: 2023-02-06
+date-format: "MMMM D, YYYY"
+description: "Bayesian stopping criteria for QMC integration and efficient credible-interval computation."
+categories:
+ - "QMCPy"
+ - "Bayesian Methods"
+ - "Stopping Criteria"
+---
+
+This post explains Bayesian stopping criteria for QMC integration and how matching kernels with lattice or digital net designs reduces credible-interval costs.
+
+The blog [Why Add Q to MC?](../why-add-q-to-mc/index.qmd) explained the advantages of carefully chosen, low discrepancy sampling sites for approximating multivariate integrals, or equivalently, expectations of functions of multivariate random variables. This blog post explains a Bayesian approach to determining the sample size required to satisfy the user's error tolerance.
+
+Recall that the problem of interest takes the following form:
+
+$$
+\begin{aligned}
+\mu
+&:= \int_{[0,1]^d} f(\boldsymbol{x}) \,
+\mathrm{d}\boldsymbol{x}
+= \mathbb{E}[f(\boldsymbol{X})],
+\qquad
+\boldsymbol{X} \sim \mathcal{U}[0,1]^d, \\
+\mu \approx \widehat{\mu}_n
+&:= \frac{1}{n}
+\sum_{i=1}^n f(\boldsymbol{x}_i),
+\qquad
+\boldsymbol{x}_1, \boldsymbol{x}_2, \ldots
+\text{ are the sampling sites or nodes.}
+\end{aligned}
+$$
+
+The question of when enough samples have been taken to satisfy the
+specified error criterion,
+$\left|\mu - \widehat{\mu}_n\right| \le \varepsilon$, is an important
+one. Stopping criteria based on Bayesian credible intervals are one good
+way of answering this question.
+
+Traditionally, Bayesian credible intervals take $\mathcal{O}(n^3)$
+operations to construct, where $n$ is the sample size. This is due to
+the need to invert an $n \times n$ covariance matrix. In contrast,
+computing $\widehat{\mu}$ requires only $\mathcal{O}(n)$ operations.
+Fortuitously, there is a way to reduce the computational cost of
+constructing the Bayesian credible interval to
+$\mathcal{O}(n \log(n))$ when using lattice or digital net low
+discrepancy samples. These Bayesian stopping criteria are implemented as
+[`CubBayesLatticeG`](https://qmcsoftware.github.io/QMCSoftware/) and
+[`CubBayesNetG`](https://qmcsoftware.github.io/QMCSoftware/). This blog
+explains the key points.
+
+## The Bayesian Approach
+
+The Bayesian approach to approximating integrals assumes that the
+integrand, $f:[0,1]^d \to \mathbb{R}$, is drawn from a Gaussian
+stochastic process parameterized by a constant mean and a covariance
+function defined by scale and shape parameters. This is denoted
+$f \sim \mathcal{GP}(m, s^2 C_{\boldsymbol{\theta}})$, and means that
+
+$$
+\mathbb{E}[f(\boldsymbol{x})] = m,
+\qquad
+\operatorname{cov}
+\left(f(\boldsymbol{t}), f(\boldsymbol{x})\right)
+=
+\mathbb{E}
+\left[
+\left(f(\boldsymbol{t}) - m\right)
+\left(f(\boldsymbol{x}) - m\right)
+\right]
+=
+s^2 C_{\boldsymbol{\theta}}(\boldsymbol{t}, \boldsymbol{x}),
+\qquad
+\forall \boldsymbol{t}, \boldsymbol{x} \in [0,1]^d.
+$$
+
+Here, $m$, $s$, and $C_{\boldsymbol{\theta}}$ must be specified or
+estimated. Because $f$ is a Gaussian process, the multivariate integral
+of $f$, i.e., $\mu$, has a Gaussian distribution. Automatic Bayesian
+cubature uses these assumptions to increase $n$ until we reach
+
+$$
+\mathbb{P}_f
+\left[
+\left| \mu - \widehat{\mu}_n \right|
+\le \varepsilon
+\right]
+\ge 99\%.
+$$
+
+The [Bayesian credible interval](https://arxiv.org/abs/1809.09803), which depends on the sampling nodes and the parameters defining the Gaussian process, is
+
+$$
+\mathbb{P}_f
+\left[
+\left|\mu - \widehat{\mu}_n\right|
+\le \operatorname{err}_{\operatorname{CI}}
+\right]
+= 99\%,
+\qquad
+\operatorname{err}_{\operatorname{CI}}
+:=
+2.58\,s
+\sqrt{
+c_{0,\boldsymbol{\theta}}
+-
+\boldsymbol{c}_{\boldsymbol{\theta}}(\mathsf{X})^{\mathsf{T}}
+\mathsf{C}_{\boldsymbol{\theta}}(\mathsf{X}, \mathsf{X})^{-1}
+\boldsymbol{c}_{\boldsymbol{\theta}}(\mathsf{X})
+},
+$$
+
+assuming that the mean of the Gaussian process, $m$, is zero. The
+formulae for the case of general $m$ are more complicated, but similar.
+Here,
+
+$$
+\begin{gathered}
+\mathsf{X}
+=
+\left(
+\boldsymbol{x}_1, \ldots, \boldsymbol{x}_n
+\right)^{\mathsf{T}},
+\qquad
+c_{0,\boldsymbol{\theta}}
+:=
+\int_{[0,1]^d \times [0,1]^d}
+C_{\boldsymbol{\theta}}(\boldsymbol{t}, \boldsymbol{x})
+\, \mathrm{d}\boldsymbol{t}
+\, \mathrm{d}\boldsymbol{x}, \\
+\boldsymbol{c}_{\boldsymbol{\theta}}(\mathsf{X})
+:=
+\left(
+\int_{[0,1]^d}
+C_{\boldsymbol{\theta}}(\boldsymbol{t}, \boldsymbol{x}_i)
+\, \mathrm{d}\boldsymbol{t}
+\right)_{i=1}^n,
+\qquad
+\mathsf{C}_{\boldsymbol{\theta}}(\mathsf{X}, \mathsf{X})
+=
+\left(
+C_{\boldsymbol{\theta}}(\boldsymbol{x}_i, \boldsymbol{x}_j)
+\right)_{i,j=1}^n,
+\end{gathered}
+$$
+
+and $2.58$ represents the $99.5\%$ quantile of the standard Gaussian
+distribution. Bayesian cubature proceeds by increasing $n$ until
+$\operatorname{err}_{\operatorname{CI}}$ is no greater than the error
+tolerance.
+
+The Bayesian credible interval as a stopping criterion assumes that the
+integrand, $f$, is not an outlier with respect to this Gaussian process.
+To improve the chances that this is the case, the hyperparameters $m$,
+$s$, and $\boldsymbol{\theta}$ should be tuned using the function data.
+In particular, $s$ should depend on the magnitude of the fluctuations of
+$f$. The standard approaches for tuning the hyperparameters require
+optimization, for which each iteration involves an eigendecomposition of
+$C_{\boldsymbol{\theta}}(\mathsf{X}, \mathsf{X})$.
+
+## Fast Bayesian Cubature
+
+The computation of $\operatorname{err}_{\operatorname{CI}}$ for a single
+$\mathsf{X}$ and $\boldsymbol{\theta}$ requires inversion of the
+covariance matrix
+$C_{\boldsymbol{\theta}}(\mathsf{X}, \mathsf{X})$, which typically costs
+$\mathcal{O}(n^3)$. An eigendecomposition of
+$C_{\boldsymbol{\theta}}(\mathsf{X}, \mathsf{X})$ costs at least this
+much. The way to overcome this computational burden is to use families
+of covariance kernels, $C_{\boldsymbol{\theta}}$, that match the design,
+$\mathsf{X}$, in a way that the eigenvectors and eigenvalues of the
+matrix $C_{\boldsymbol{\theta}}(\mathsf{X}, \mathsf{X})$ can be obtained
+via fast transforms.
+
+For
+[integration lattices](https://global.oup.com/academic/product/lattice-methods-for-multiple-integration-9780198534723?cc=us&lang=en&)
+and kernels $C_{\boldsymbol{\theta}}$ such as
+
+$$
+C_{\boldsymbol{\theta}}(\boldsymbol{t}, \boldsymbol{x})
+=
+K_{\boldsymbol{\theta}}
+\left(
+\boldsymbol{t} - \boldsymbol{x} \bmod \boldsymbol{1}
+\right),
+$$
+
+where
+
+$$
+K_{\boldsymbol{\theta}}(\boldsymbol{x})
+=
+\prod_{l=1}^d
+\left[
+1
+- (-1)^r \eta B_{2r}(x_l)
+\right],
+\qquad
+\boldsymbol{\theta} = (r,\eta),
+\quad
+r \in \mathbb{N},
+\quad
+\eta > 0,
+$$
+
+and $B_{2r}$ is an even degree Bernoulli polynomial, the fast transform
+corresponds to a fast
+[Fourier transform](https://www.ams.org/journals/mcom/1965-19-090/S0025-5718-1965-0178586-1/S0025-5718-1965-0178586-1.pdf).
+For digital nets, see the blog
+[Digital Sequences, the Niederreiter Construction](https://qmcpy.org/2021/06/04/digital-sequences-the-niederreiter-construction/).
+For $C_{\boldsymbol{\theta}}$ that are
+[digitally-shift invariant](https://repository.iit.edu/islandora/object/islandora%3A1009768),
+the fast transform corresponds to a
+[fast Walsh transform](https://dl.acm.org/doi/abs/10.1109/TC.1976.1674569).
+In both cases, the computational burden attributable to tuning the
+hyperparameters and computing the width of the credible interval is a
+reasonable $\mathcal{O}(n \log(n))$.
+
+The following code shows how to use a Bayesian stopping criterion for Keister's example.
+
+```python
+import qmcpy as qp
+
+tol = 0.005
+
+integrand = qp.Keister(qp.Lattice(dimension=2, order="RADICAL INVERSE"))
+keister_2d_exact = integrand.exact_integ(2)
+solution, data = qp.CubBayesLatticeG(
+ integrand,
+ abs_tol=tol,
+ n_init=2**5,
+).integrate()
+print("Integration error: ", abs(solution - keister_2d_exact) < tol)
+```
+
+Listing 1: Example usage of the lattice Bayesian cubature algorithm.
+
+This example can be run in Google Colab without any installation using this [notebook](https://github.com/QMCSoftware/QMCSoftware/blob/develop/demos/integration_examples.ipynb).
+
+## References
+
+1. Hickernell, F. J. Blog: Why Add Q to MC?
+[https://qmcpy.org/2020/06/25/why_add_q_to_mc/](https://qmcpy.org/2020/06/25/why_add_q_to_mc/). 2020.
+2. Choi, S.-C. T., Hickernell, F., McCourt, M., & Sorokin, A. QMCPy: A
+quasi-Monte Carlo Python Library. [https://qmcsoftware.github.io/QMCSoftware/](https://qmcsoftware.github.io/QMCSoftware/). 2020.
+3. Rathinavel, J., & Hickernell, F. Fast automatic Bayesian cubature
+using lattice sampling. *Statistics and Computing*, 29, 1215-1229 (2019).
+4. Sloan, I. H., & Joe, S. *Lattice Methods for Multiple Integration*.
+ Oxford University Press, Oxford (1994).
+5. Hickernell, F., & Niederreiter, H. The existence of good extensible
+ rank-1 lattices. *Journal of Complexity*, 19, 286-300 (2003).
+6. Cooley, J. W., & Tukey, J. W. An algorithm for the machine
+calculation of complex Fourier series. *Mathematics of Computation*, 19, 297-301 (1965).
+7. Ebert, A. Blog: Digital Sequences, the Niederreiter Construction.
+[https://qmcpy.org/2021/06/04/digital-sequences-the-niederreiter-construction/](https://qmcpy.org/2021/06/04/digital-sequences-the-niederreiter-construction/). 2021.
+8. Rathinavel, J. *Fast Automatic Bayesian Cubature Using Matching
+Kernels and Designs*. PhD thesis, Illinois Institute of Technology (2019).
+9. Fino, B. J., & Algazi, V. R. Unified matrix treatment of the fast
+Walsh-Hadamard transform. *IEEE Transactions on Computers*, C-25, 1142-1146 (1976).
diff --git a/blog/cubmccltvec-vectorizing-the-cubmcclt-algorithm/index.qmd b/blog/cubmccltvec-vectorizing-the-cubmcclt-algorithm/index.qmd
new file mode 100644
index 0000000..0f88d27
--- /dev/null
+++ b/blog/cubmccltvec-vectorizing-the-cubmcclt-algorithm/index.qmd
@@ -0,0 +1,232 @@
+---
+title: "CubMCCLTVec: Vectorizing the CubMCCLT Algorithm"
+author: "Aadit Jain"
+date: 2026-02-25
+date-format: "MMMM D, YYYY"
+description: "A vectorized extension of CubMCCLT for confidence intervals on vector-valued quantities of interest."
+categories:
+ - "QMCPy"
+ - "Vector-Valued Integration"
+ - "Stopping Criteria"
+---
+
+This post introduces `CubMCCLTVec`, a vectorized extension of `CubMCCLT` for confidence intervals on vector-valued quantities of interest.
+
+Recent work by Aleksei G. Sorokin and Jagadeeswaran Rathinavel [1]
+discuss extending stopping criterion for a scalar mean to stopping
+criterion for vector quantities of interest formulated as functions of
+multiple means. One such scalar stopping criterion is `CubMCCLT` that
+calculates a confidence interval for $\mu$ by using the Central Limit
+Theorem. When $\{x_i\}_{i=1}^{n}$ are IID and $f$ has a finite
+variance, the confidence interval for $\mu$ can be calculated as:
+
+$$
+\mu^{\pm} = \hat{\mu}\; \pm\; C Z_{\alpha^{(\mu)}/2} \hat{\sigma}/\sqrt{n}
+$$
+
+Here, $\hat{\mu}$ is the sample average of function evaluations,
+$Z_{\alpha^{(\mu)}/2}$ is the inverse CDF of a standard normal
+distribution at $1 - \alpha^{(\mu)}/2$, the variance of $f(X)$ can
+be approximated by
+$\hat{\sigma}^2 = \frac{1}{n-1} \sum_{i=1}^{n} (f(\boldsymbol{x}_i) - \hat{\mu})^2$,
+and $C^2$ is an inflation factor greater than 1 for a more
+conservative estimate.
+
+Building on the `CubMCCLT` algorithm, we have developed a vectorized version of it known as `CubMCCLTVec`.
+
+## What Does the CubMCCLTVec Class Do?
+
+The `CubMCCLTVec` class, which is a stopping criterion object, calculates a confidence interval for functions with multiple outputs based on the user-specified confidence level (default is 99%). Given an initial and maximum sample size and an absolute tolerance, we keep on doubling the sample size and recomputing the confidence interval until half the confidence interval width is less than the absolute tolerance or the double of the current sample size exceeds the maximum sample size.
+
+Like the other stopping criterion objects, `CubMCCLTVec` too utilizes an accumulate data object to recompute the confidence interval known as `MeanVarDataVec`.
+
+## Some Examples Utilizing the CubMCCLTVec Class
+
+The following code illustrates three examples that are being solved using `CubMCCLTVec`:
+
+1. Covariance [2]: $T \sim \mathcal{N}(1,I_d)$ and the covariance of
+ $P = T_0\cdots T_{d-1}$ and $S = T_0+\dots+T_{d-1}$ is:
+
+ $$
+ \mathrm{Cov}[P,S] = \mathbb{E}[PS]-\mathbb{E}[P]\mathbb{E}[S] = \mu_0-\mu_1\mu_2
+ $$
+
+ Theoretically, $\mathrm{Cov}[P,S] = 2d-(1)(d) = d$.
+
+2. Box Integral [3]: $B_s(x) = (\sum_{j=1}^{d} x^2_j)^{s/2}$ where
+ $x_1,\dots,x_d \sim \mathcal{U}[0,1]$ and the box integral is
+ computed for $s=-1$ and $s=1$ (the two outputs).
+
+3. Custom Fun: $x_j \sim \mathcal{U}[0,2j]$ for $j=1,\dots,6$.
+
+## Python Implementation
+
+```python
+import qmcpy as qp
+import numpy as np
+
+# Example 1: Covariance
+dimension = 4
+true_measure = qp.IIDStdUniform(dimension)
+
+class Covariance(qp.integrand.Integrand):
+ def __init__(self, true_measure):
+ super().__init__(true_measure)
+
+ def g(self, x):
+ P = np.prod(x, axis=1)
+ S = np.sum(x, axis=1)
+ return np.vstack((P*S, P, S)).T
+
+integrand = Covariance(true_measure)
+solution = qp.CubMCCLTVec(integrand, abs_tol=1e-2)
+print("Covariance:", solution.integrate())
+
+
+# Example 2: Box Integral
+dimension = 5
+true_measure = qp.IIDStdUniform(dimension)
+
+class BoxIntegral(qp.integrand.Integrand):
+ def __init__(self, true_measure):
+ super().__init__(true_measure)
+
+ def g(self, x):
+ norm_sq = np.sum(x**2, axis=1)
+ return np.vstack((norm_sq**(-0.5), norm_sq**(0.5))).T
+
+integrand = BoxIntegral(true_measure)
+solution = qp.CubMCCLTVec(integrand, abs_tol=1e-2)
+print("Box Integral:", solution.integrate())
+
+
+# Example 3: Custom Function
+dimension = 6
+true_measure = qp.IIDStdUniform(dimension)
+
+class CustomFun(qp.integrand.Integrand):
+ def __init__(self, true_measure):
+ super().__init__(true_measure)
+
+ def g(self, x):
+ weights = np.arange(1, dimension+1)
+ scaled = x * (2*weights)
+ f1 = np.sum(scaled, axis=1)
+ f2 = np.prod(scaled, axis=1)
+ return np.vstack((f1, f2)).T
+
+integrand = CustomFun(true_measure)
+solution = qp.CubMCCLTVec(integrand, abs_tol=1e-2)
+print("Custom Fun:", solution.integrate())
+```
+
+**Example 1 Output: Covariance**
+
+```text
+Solution: 3.998147791818197
+Data:
+MeanVarDataVec (AccumulateData Object)
+ solution 3.998
+ comb_bound_low 3.977
+ comb_bound_high 4.019
+ comb_flags 1
+ n_total 2^(26)
+ n [67108864. 67108864. 67108864.]
+ time_integrate 152.068
+CubMCCLTVec (StoppingCriterion Object)
+ inflate 1.200
+ alpha 0.010
+ abs_tol 0.025
+ rel_tol 0
+ n_init 2^(8)
+ n_max 2^(30)
+CovIntegrand (Integrand Object)
+Gaussian (TrueMeasure Object)
+ mean 1
+ covariance 1
+ decomp_type PCA
+IIDStdUniform (DiscreteDistribution Object)
+ d 2^(2)
+ entropy 7
+ spawn_key ()
+```
+
+**Example 2 Output: Box Integral**
+
+```text
+Solution: [1.1853359 0.95670595]
+Data:
+MeanVarDataVec (AccumulateData Object)
+ solution [1.185 0.957]
+ comb_bound_low [1.139 0.918]
+ comb_bound_high [1.232 0.995]
+ comb_flags [ True True]
+ n_total 2^(11)
+ n [2048. 512.]
+ time_integrate 0
+CubMCCLTVec (StoppingCriterion Object)
+ inflate 1.200
+ alpha 0.010
+ abs_tol 0.050
+ rel_tol 0
+ n_init 2^(8)
+ n_max 2^(30)
+BoxIntegral (Integrand Object)
+ s [-1 1]
+Uniform (TrueMeasure Object)
+ lower_bound 0
+ upper_bound 1
+IIDStdUniform (DiscreteDistribution Object)
+ d 3
+ entropy 7
+ spawn_key ()
+```
+
+**Example 3 Output: Custom Fun**
+
+```text
+Solution: [[1.00001237 1.99884832 2.99870039]
+ [4.00056276 4.99789845 6.00074045]]
+Data:
+MeanVarDataVec (AccumulateData Object)
+ solution [[1. 1.999 2.999]
+ [4.001 4.998 6.001]]
+ comb_bound_low [[0.99 1.989 2.991]
+ [3.991 4.989 5.993]]
+ comb_bound_high [[1.01 2.009 3.006]
+ [4.01 5.007 6.008]]
+ comb_flags [[ True True True]
+ [ True True True]]
+ n_total 2^(21)
+ n [[ 32768. 131072. 524288.]
+ [ 524288. 1048576. 2097152.]]
+ time_integrate 0.650
+CubMCCLTVec (StoppingCriterion Object)
+ inflate 1.200
+ alpha 0.010
+ abs_tol 0.010
+ rel_tol 0
+ n_init 2^(8)
+ n_max 2^(30)
+CustomFun (Integrand Object)
+Uniform (TrueMeasure Object)
+ lower_bound 0
+ upper_bound 1
+IIDStdUniform (DiscreteDistribution Object)
+ d 6
+ entropy 7
+ spawn_key ()
+```
+
+## Benefits of Developing the CubMCCLTVec Class
+
+This class gives us a new and different option to find when the user-specified error tolerance has been satisfied and its generalization to functions with multiple outputs allows us to utilize the existing `CubMCCLT` algorithm and extend it to such functions.
+
+## References
+
+1. Sorokin, A. G., & Rathinavel, J. *On Bounding and Approximating
+Functions of Multiple Expectations using Quasi-Monte Carlo*. To appear in the *Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing Proceedings 2022* (2022).
+2. Sorokin, A. *Monte Carlo for Vector Functions of Integrals*. Jupyter
+Notebook. QMCPy: A quasi-Monte Carlo Python Library, 2023. [https://github.com/QMCSoftware/QMCSoftware/blob/master/demos/pydata.chi.2023.ipynb](https://github.com/QMCSoftware/QMCSoftware/blob/master/demos/pydata.chi.2023.ipynb).
+3. Bailey, D., Borwein, J., & Crandall, R. *Box integrals*. *Journal of
+Computational and Applied Mathematics* **206**, 196-208. ISSN: 0377-0427. [https://www.sciencedirect.com/science/article/pii/S0377042706004250](https://www.sciencedirect.com/science/article/pii/S0377042706004250) (2007).
diff --git a/blog/digital-sequences-the-niederreiter-construction/figures/cubqmcnetg-samples.png b/blog/digital-sequences-the-niederreiter-construction/figures/cubqmcnetg-samples.png
new file mode 100644
index 0000000..5549ebe
Binary files /dev/null and b/blog/digital-sequences-the-niederreiter-construction/figures/cubqmcnetg-samples.png differ
diff --git a/blog/digital-sequences-the-niederreiter-construction/figures/keister-error.png b/blog/digital-sequences-the-niederreiter-construction/figures/keister-error.png
new file mode 100644
index 0000000..f681b1f
Binary files /dev/null and b/blog/digital-sequences-the-niederreiter-construction/figures/keister-error.png differ
diff --git a/blog/digital-sequences-the-niederreiter-construction/figures/niederreiter-points.png b/blog/digital-sequences-the-niederreiter-construction/figures/niederreiter-points.png
new file mode 100644
index 0000000..32a488c
Binary files /dev/null and b/blog/digital-sequences-the-niederreiter-construction/figures/niederreiter-points.png differ
diff --git a/blog/digital-sequences-the-niederreiter-construction/index.qmd b/blog/digital-sequences-the-niederreiter-construction/index.qmd
new file mode 100644
index 0000000..9aa1907
--- /dev/null
+++ b/blog/digital-sequences-the-niederreiter-construction/index.qmd
@@ -0,0 +1,301 @@
+---
+title: "Digital Sequences, the Niederreiter Construction"
+author: "Adrian Ebert"
+date: 2021-06-07
+date-format: "MMMM D, YYYY"
+description: "The digital construction of Niederreiter sequences and their performance relative to Sobol' sequences."
+categories:
+ - "QMCPy"
+ - "Low-Discrepancy Sequences"
+ - "Quasi-Monte Carlo"
+image: figures/niederreiter-points.png
+---
+
+This post explains the digital construction of Niederreiter sequences and compares their QMCPy performance with Sobol' sequences.
+
+The previous blog post on
+[What Makes a Sequence "Low Discrepancy"?](../what-makes-a-sequence-low-discrepancy/index.qmd)
+introduced the concept of so-called low discrepancy (LD) points. In the
+literature on QMC methods, there are in general two main families of low
+discrepancy point sets that are commonly used as integration nodes.
+These are, on the one hand, lattice point sets, as introduced
+independently by Korobov and Hlawka, and, on the other hand, digital
+$(t,m,d)$-nets and $(t,d)$-sequences, as introduced by Niederreiter,
+building up on ideas by Sobol' and Faure. In this post, we will take a
+closer look at digital nets and sequences. In particular, we will
+investigate the so-called *Niederreiter sequence*, as introduced in [4].
+
+## Digital Construction Scheme
+
+We introduce the construction scheme of digital $(t,d)$-nets and
+sequences. Here, the dimension of the point set is denoted by $d$ and
+$t$ is the so-called quality parameter of the digital sequence or net.
+The resulting point sets in the unit cube $[0,1]^d$ can, for example,
+be used for quasi-Monte Carlo integration. For a more detailed overview
+of digital nets and sequences, we refer to the survey article [1] or the
+monograph [2].
+
+We will use the following notation. Let $b$ be prime and let
+$\mathbb{Z}_b := \{\overline{0}, \overline{1}, \dots,
+\overline{b-1}\}$ denote the equivalence classes of integers modulo
+$b$. Equipped with addition and multiplication, $\mathbb{Z}_b$ is a
+finite field of order $b$. With a slight abuse of notation, we will
+identify the elements of the field $\mathbb{Z}_b$ with the integers
+$\{0,1,\ldots,b-1\}$.
+
+Given the matrices
+$C_1,\dots,C_d \in \mathbb{Z}_b^{\mathbb{N}\times\mathbb{N}}$ with
+$C_j = (c_{j,k,\ell})_{k,\ell \in \mathbb{N}}$, the $j$-th component of
+the $i$-th point of the sequence
+$S = (\boldsymbol{t}_0,\boldsymbol{t}_1,\dots) \in [0,1]^d$, denoted by
+$t_{i,j}$, is then constructed as follows:
+
+1. Write $i \in \mathbb{N}$ in its base $b$ representation, that is,
+
+ $$
+ i = (\cdots i_3 i_2 i_1)_b
+ = i_1 + i_2 b + i_3 b^2 + \cdots.
+ $$
+
+2. Compute the matrix-vector product
+
+ $$
+ \begin{pmatrix}
+ y_1\\
+ y_2\\
+ y_3\\
+ \vdots
+ \end{pmatrix}
+ =
+ C_j
+ \begin{pmatrix}
+ i_1\\
+ i_2\\
+ i_3\\
+ \vdots
+ \end{pmatrix},
+ $$
+
+ where all additions and multiplications are done over the finite
+ field $\mathbb{Z}_b$.
+
+3. Set the $j$-th component of the $i$-th point as
+
+ $$
+ t_{i,j}
+ = (0.y_1 y_2 \dots)_b
+ = \frac{y_1}{b} + \frac{y_2}{b^2}
+ + \frac{y_3}{b^3} + \cdots.
+ $$
+
+The resulting sequence $S = (\boldsymbol{t}_0,\boldsymbol{t}_1,\dots)$
+is then called a *digital sequence* over the finite field
+$\mathbb{Z}_b$, and the matrices $C_1,\dots,C_d$ are called the
+*generating matrices* of the digital sequence $S$.
+
+## The Niederreiter Construction
+
+A special instance of digital sequences is the so-called Niederreiter
+sequence, which was the first known construction yielding
+$(t,d)$-sequences for arbitrary dimensions $d$ and arbitrary prime
+powers $b$. The construction method is based on polynomial arithmetic
+over finite fields as we will illustrate below. Further details can be
+found in [5] and [2].
+
+Denote the set of polynomials over the finite field $\mathbb{Z}_b$ by
+$\mathbb{Z}_b[x]$. Furthermore, let $\mathbb{Z}_b((x^{-1}))$ be the
+field of formal Laurent series over $\mathbb{Z}_b$. Elements of
+$\mathbb{Z}_b((x^{-1}))$ are of the form
+
+$$
+L = \sum_{\ell = w}^{\infty} t_\ell \, x^{-\ell},
+$$
+
+where $w$ is an arbitrary integer and the coefficients $t_\ell$ are
+elements of $\mathbb{Z}_b$. The set of formal Laurent series contains
+the set $\mathbb{Z}_b[x]$ as a subset and also permits addition,
+subtraction, multiplication, and division of formal series.
+
+We summarize the construction scheme for the generating matrices of the
+Niederreiter sequence, which is based on arithmetic in the field of
+Laurent series $\mathbb{Z}_b((x^{-1}))$, below.
+
+1. Let $p_1,\ldots,p_d \in \mathbb{Z}_b[x]$ be distinct monic
+ irreducible polynomials over $\mathbb{Z}_b$. For $1 \le j \le d$,
+ let $e_j := \deg(p_j)$ be the degree of the polynomial $p_j$.
+
+2. For integers $1 \le j \le d$, $r \ge 1$, and $0 \le k < e_j$,
+ consider the expansions
+
+ $$
+ \frac{x^k}{p_j(x)^r}
+ =
+ \sum_{\ell = 0}^{\infty}
+ a^{(j)}(r,k,\ell) \, x^{-\ell-1}
+ $$
+
+ over the field of formal Laurent series
+ $\mathbb{Z}_b((x^{-1}))$.
+
+3. We then define the entries of the matrix
+ $C_j = (c_{j,i,\ell})_{i \ge 1,\ell \ge 0}$ as
+
+ $$
+ c_{j,i,\ell}
+ :=
+ a^{(j)}(Q+1,k,\ell)
+ \in \mathbb{Z}_b
+ \qquad
+ \text{for } 1 \le j \le d,\; i \ge 1,\; \ell \ge 0,
+ $$
+
+ where $i-1 = Q e_j + k$ with integers $Q=Q(i,j)$ and $k=k(i,j)$
+ with $0 \le k < e_j$.
+
+The digital sequence which is generated by the resulting generating
+matrices $C_1,\ldots,C_d$ is called the Niederreiter sequence. In
+particular, it is a digital $(t,d)$-sequence over the finite field
+$\mathbb{Z}_b$ and the quality parameter $t$ equals
+
+$$
+t = \sum_{j=1}^d (e_j - 1).
+$$
+
+The $t$-value of a digital net or sequence is a quality parameter of
+uniformity for the corresponding point set. A digital sequence is well
+distributed if the $t$-value is small.
+
+In order to illustrate the structure of the Niederreiter points, we
+display the first $128$ points of the sequence in 2 dimensions in Figure
+1 below.
+
+
+
+ Figure 1: The first $128$ points of the Niederreiter sequence in 2 dimensions.
+
+
+The performance of the Niederreiter sequence in practical applications is, in general, similar to that of the widely used Sobol' sequence. For a more detailed comparison in financial applications, see for example [3].
+
+## Niederreiter Points via QMCPy
+
+The Niederreiter sequence has recently been added to the `DiscreteDistribution` class of the QMCPy Python library. As a digital sequence, the Niederreiter sequence is part of the `DigitalNet` or `Sobol` generator and can be accessed by specifying the corresponding generating matrices. The code in Listing 1 below can be used to draw randomized points from QMCPy's Niederreiter object.
+
+```python
+from qmcpy import DigitalNet
+
+nied = DigitalNet(dimension=5, z_path="niederreiter_mat.20000.32.msb.npy")
+x = nied.gen_samples(n_min=0, n_max=4)
+x
+```
+
+```python
+array([[0.89 , 0.603, 0.288, 0.881, 0.298],
+ [0.224, 0.063, 0.052, 0.854, 0.387],
+ [0.707, 0.266, 0.743, 0.739, 0.07 ],
+ [0.43 , 0.806, 0.98 , 0.527, 0.246]])
+```
+
+Listing 1: Supplying the Niederreiter generating matrices to the digital net generator.
+
+For further details on the construction of the generating matrices of the Niederreiter sequence, we refer the interested reader to the [documentation](https://bitbucket.org/adrian_ebert/qmc-construction-algorithms/src/master/digital-constructions/digseq/niederreitermats/construction_programmes/documentation-niederreiter.pdf). Additionally, the C++ code which was used to construct the generating matrices can be found [here](https://bitbucket.org/adrian_ebert/qmc-construction-algorithms/src/master/digital-constructions/digseq/niederreitermats/construction_programmes/niederreiter.cpp).
+
+## Comparison with the Sobol' Sequence
+
+In order to demonstrate the effectiveness of the Niederreiter sequence
+for numerical integration, we consider its performance for a test
+function and compare it with the commonly used Sobol' sequence. In
+particular, we consider the problem of integrating the Keister function
+with respect to a $d$-dimensional Gaussian measure:
+
+$$
+\begin{aligned}
+f(\boldsymbol{x})
+&=
+\pi^{d/2}
+\cos\left(\|\boldsymbol{x}\|\right),
+\qquad
+\boldsymbol{x} \in \mathbb{R}^d,
+\qquad
+\boldsymbol{X} \sim
+\mathcal{N}(\boldsymbol{0}, \boldsymbol{I}_d/2), \\
+I(f)
+&=
+\int_{[0,1]^d}
+\pi^{d/2}
+\cos\left(
+\sqrt{
+\frac{1}{2}
+\sum_{j=1}^d
+\left(\Phi^{-1}(x_j)\right)^2
+}
+\right)
+\, \mathrm{d}\boldsymbol{x},
+\end{aligned}
+$$
+
+where $\|\cdot\|$ is the Euclidean norm, $\boldsymbol{I}_d$ is the
+$d$-dimensional identity matrix, and $\Phi$ denotes the standard normal
+cumulative distribution function.
+
+The value of $I(f)$ is then approximated by a QMC rule which uses
+randomized Sobol' or Niederreiter points as cubature nodes, i.e.,
+
+$$
+I(g)
+\approx
+Q_N(g)
+=
+\frac{1}{N}
+\sum_{n=0}^{N-1}
+g\left(\boldsymbol{x}_n\right),
+$$
+
+where
+
+$$
+g(\boldsymbol{x})
+=
+\pi^{d/2}
+\cos\left(
+\sqrt{
+\frac{1}{2}
+\sum_{j=1}^d
+\left(\Phi^{-1}(x_j)\right)^2
+}
+\right).
+$$
+
+In Figure 2 we display the mean-square error of the approximation
+$Q_N(g)$ for different numbers $N$ of total function evaluations and a
+dimension of $d=5$. Here we have computed the exact solution and track
+the decay of the error as we increase the sample size. For the Keister
+test function, the numerical results show that the Niederreiter sequence
+is competitive with the Sobol' sequence as both errors decay at an
+almost identical rate.
+
+
+
+ Figure 2: Error convergence behavior for the numerical integration of the Keister function with $d=5$ using Niederreiter and Sobol' points.
+
+
+QMCPy also includes stopping criteria that automatically select the number of points required to meet an error tolerance. Here it is assumed that the exact solution is not known and we track the number of samples required to guarantee an approximation within a user-specified tolerance. Figure 3 shows similar performance for Niederreiter and Sobol' sequences when utilized by the `CubQMCNetG` stopping criterion.
+
+
+
+ Figure 3: Number of samples required to meet user-specified error tolerances using the CubQMCNetG stopping criterion with Niederreiter and Sobol' sequences for the $d=5$ Keister function.
+
+
+## References
+
+1. Dick, J., Kuo, F. Y., & Sloan, I. H. High-dimensional integration:
+ The quasi-Monte Carlo way. *Acta Numerica*, 22, 133-288 (2013).
+2. Dick, J., & Pillichshammer, F. *Digital Nets and Sequences*.
+ Cambridge University Press (2010).
+3. Harase, S. Comparison of Sobol' sequences in financial applications.
+ *Monte Carlo Methods and Applications*, 25(1), 61-74 (2019).
+4. Niederreiter, H. Point sets and sequences with small discrepancy.
+ *Monatshefte fur Mathematik*, 104, 273-337 (1987).
+5. Niederreiter, H. *Random Number Generation and Quasi-Monte Carlo
+Methods*. Number 63 in CBMS-NSF Series in Applied Mathematics. SIAM, Philadelphia (1992).
+6. Bitbucket repository of Adrian Ebert.
+ [https://bitbucket.org/adrian_ebert/qmc-construction-algorithms/](https://bitbucket.org/adrian_ebert/qmc-construction-algorithms/).
diff --git a/blog/gbm-qmcpy/figures/figure_1.png b/blog/gbm-qmcpy/figures/figure_1.png
new file mode 100644
index 0000000..b15df11
Binary files /dev/null and b/blog/gbm-qmcpy/figures/figure_1.png differ
diff --git a/blog/gbm-qmcpy/figures/figure_2.png b/blog/gbm-qmcpy/figures/figure_2.png
new file mode 100644
index 0000000..1073cde
Binary files /dev/null and b/blog/gbm-qmcpy/figures/figure_2.png differ
diff --git a/blog/gbm-qmcpy/figures/figure_3.png b/blog/gbm-qmcpy/figures/figure_3.png
new file mode 100644
index 0000000..f62c638
Binary files /dev/null and b/blog/gbm-qmcpy/figures/figure_3.png differ
diff --git a/blog/gbm-qmcpy/figures/figure_4.png b/blog/gbm-qmcpy/figures/figure_4.png
new file mode 100644
index 0000000..07ae4b4
Binary files /dev/null and b/blog/gbm-qmcpy/figures/figure_4.png differ
diff --git a/blog/gbm-qmcpy/figures/figure_5.png b/blog/gbm-qmcpy/figures/figure_5.png
new file mode 100644
index 0000000..574a254
Binary files /dev/null and b/blog/gbm-qmcpy/figures/figure_5.png differ
diff --git a/blog/gbm-qmcpy/figures/figure_6.png b/blog/gbm-qmcpy/figures/figure_6.png
new file mode 100644
index 0000000..3dfdbb2
Binary files /dev/null and b/blog/gbm-qmcpy/figures/figure_6.png differ
diff --git a/blog/gbm-qmcpy/figures/figure_7.png b/blog/gbm-qmcpy/figures/figure_7.png
new file mode 100644
index 0000000..67f44f6
Binary files /dev/null and b/blog/gbm-qmcpy/figures/figure_7.png differ
diff --git a/blog/gbm-qmcpy/index.qmd b/blog/gbm-qmcpy/index.qmd
new file mode 100644
index 0000000..be35cf0
--- /dev/null
+++ b/blog/gbm-qmcpy/index.qmd
@@ -0,0 +1,414 @@
+---
+title: "Highly Efficient Geometric Brownian Motion Modeling with QMCPy"
+author: "Larysa Matiukha, Aleksei Sorokin, and Sou-Cheng Choi"
+date: 2025-10-22
+date-format: "MMMM D, YYYY"
+description: "Generating and analyzing Geometric Brownian Motion paths efficiently with QMCPy."
+categories:
+ - "QMCPy"
+ - "Finance"
+ - "Geometric Brownian Motion"
+image: figures/figure_2.png
+---
+
+This post demonstrates how QMCPy can generate and analyze Geometric Brownian Motion paths efficiently for simulation and finance workflows.
+
+::: {.callout-note}
+## Example implementation of GBM using QMCPy
+
+Install the required Python packages:
+
+```bash
+pip install qmcpy numpy matplotlib
+```
+
+Generate 16 paths on $[0,1]$ with QMCPy's sampler ($S_0=1$,
+$\mu=0.05$, $\sigma^2=0.2$) and plot five:
+
+```python
+import numpy as np
+import matplotlib.pyplot as plt
+import qmcpy as qp
+
+sampler = qp.Lattice(252, seed=42) # daily steps for 1 year
+gbm = qp.GeometricBrownianMotion(
+ sampler, t_final=1, initial_value=1, drift=0.05, diffusion=0.2
+)
+paths = gbm.gen_samples(16)
+t = np.linspace(0, 1.0, paths.shape[1])
+plt.plot(t, paths[:5].T, alpha=0.8)
+plt.xlabel("$t$")
+plt.ylabel("$S_t$")
+plt.title("GBM paths")
+plt.show()
+```
+:::
+
+## Introduction
+
+In this blog, we demonstrate how to simulate and analyze a geometric Brownian motion (GBM) process using QMCPy in Python. GBM is widely used in finance to model stock prices and other assets. We will walk through key code snippets, plots, and insights. The numerical results can be reproduced using the [GBM demo notebook](https://qmcsoftware.github.io/QMCSoftware/demos/GBM/gbm_demo/).
+
+GBM is a continuous stochastic process in which the natural logarithm of its values follows a Brownian motion (BM) [1]. Mathematically, it can be defined as follows:
+
+$$
+S_t = S_0 \, e^{\bigl(\mu - \tfrac{\sigma^2}{2}\bigr)t + \sigma W_t}.
+$$
+
+where
+
+- $S_0$ is the initial value,
+- $\mu$ is a drift coefficient,
+- $\sigma$ is the volatility (note: QMCPy uses diffusion $= \sigma^2$),
+- $W_t$ is a standard BM.
+
+At any time $0 < t \le T$, where $0$ and $T$ represent the beginning and
+end time of the process, $S_t$ follows a log-normal distribution with
+expected value and variance as follows (see Section 3.2 in [1]):
+
+- $E[S_t] = S_0 e^{\mu t}$,
+- $\operatorname{Var}[S_t] = S_0^2 e^{2\mu t}(e^{\sigma^2 t} - 1)$,
+- $\operatorname{Cov}(S_{t_i}, S_{t_j}) = S_0^2 e^{\mu(t_i + t_j)}
+ \left(e^{\sigma^2 \min(t_i, t_j)} - 1\right)$.
+
+GBM is commonly used to model stock prices driving option payoffs in derivatives pricing [2, 3].
+
+## GBM Objects in QMCPy
+
+GBM in QMCPy inherits from `BrownianMotion` [4, 5]. We can instantiate a GBM class and generate sample paths to see the class in action:
+
+```python
+import qmcpy as qp
+
+qp_gbm = qp.GeometricBrownianMotion(qp.Lattice(2, seed=42))
+qp_gbm.gen_samples(n=4)
+```
+
+The output shows 4 sample paths evaluated at 2 time points, yielding a
+$(4 \times 2)$ array where rows represent paths and columns represent
+time steps:
+
+$$
+\left[
+\begin{array}{rr}
+0.72608046 & 0.70071241 \\
+0.38739775 & 0.07432173 \\
+0.81262942 & 1.66867239 \\
+0.61937100 & 0.31898397
+\end{array}
+\right].
+$$
+
+## Log-Normality Property
+
+The log-normal property is fundamental in financial modeling because it ensures asset prices remain strictly positive while allowing for unlimited upside potential. This property makes GBM the cornerstone of the Black-Scholes model and many derivative pricing frameworks.
+
+To validate theoretical properties, we generate $2^{12} = 4096$ paths
+over 5 time steps and compare empirical moments with theoretical values.
+The theoretical values match the last values captured in
+`qp_gbm.mean_gbm` and `qp_gbm.covariance_gbm` for the final time point.
+
+```python
+# Generate GBM samples for theoretical validation
+import qmcpy as qp
+import numpy as np
+
+S0, mu, sigma, T, n_samples = 100.0, 0.05, 0.20, 1.0, 2**12
+diffusion = sigma**2
+sampler = qp.Lattice(5, seed=42)
+qp_gbm = qp.GeometricBrownianMotion(
+ sampler, t_final=T, initial_value=S0, drift=mu, diffusion=diffusion
+)
+paths = qp_gbm.gen_samples(n_samples)
+S_T = paths[:, -1] # Final values only
+
+# Calculate theoretical vs empirical sample moments
+theo_mean = S0 * np.exp(mu * T)
+theo_var = S0**2 * np.exp(2 * mu * T) * (np.exp(diffusion * T) - 1)
+qp_emp_mean = np.mean(S_T)
+qp_emp_var = np.var(S_T, ddof=1)
+print(f"Mean: {qp_emp_mean:.3f} (theoretical: {theo_mean:.3f})")
+print(f"Variance: {qp_emp_var:.3f} (theoretical: {theo_var:.3f})")
+qp_gbm
+```
+
+| Statistic | Value |
+| --- | --- |
+| Sample Mean | 105.127 (Theoretical: 105.127) |
+| Sample Variance | 449.776 (Theoretical: 451.029) |
+| Time Vector | $[0.2,\; 0.4,\; 0.6,\; 0.8,\; 1.0]$ |
+| Drift ($\mu$) | 0.05 |
+| Diffusion ($\sigma^2$) | 0.040 |
+| Mean | $[101.005,\; 102.020,\; 103.045,\; 104.081,\; 105.127]$ |
+| Decomposition Type | PCA |
+| Covariance Matrix | $\left[\begin{array}{rrrrr} 81.943 & 82.767 & 83.599 & 84.439 & 85.288 \\ 82.767 & 167.869 & 169.556 & 171.260 & 172.981 \\ 83.599 & 169.556 & 257.923 & 260.516 & 263.134 \\ 84.439 & 171.260 & 260.516 & 352.258 & 355.798 \\ 85.288 & 172.981 & 263.134 & 355.798 & 451.029 \end{array}\right]$ |
+
+## GBM vs BM
+
+Below we compare BM and GBM using the same parameters:
+`drift = 0`, `diffusion = 1`, and `initial_value = 1`. The driftless BM
+paths should fluctuate symmetrically around the initial value ($y = 1$)
+and can take negative values, while the GBM paths remain strictly
+positive.
+
+```python
+n = 16
+sampler = qp.Lattice(2**7, seed=42)
+plot_paths("BM", sampler, t_final=1, initial_value=1, drift=0, diffusion=1, n=n)
+plot_paths("GBM", sampler, t_final=1, initial_value=1, drift=0, diffusion=1, n=n)
+```
+
+Next, we demonstrate how easily one can swap samplers or change
+parameters in QMCPy. For example, to model a stock price with initial
+value $S_0=50$, drift $\mu=0.1$, and volatility $\sigma=\sqrt{0.2}$
+over a 5-year horizon using IID sampling:
+
+```python
+gbm_iid = plot_paths(
+ "GBM",
+ qp.IIDStdUniform(2**8, seed=42),
+ t_final=5,
+ initial_value=50,
+ drift=0.1,
+ diffusion=0.2,
+ n=32,
+)
+```
+
+We can also use a low-discrepancy lattice sampler with the same parameters:
+
+```python
+gbm_lattice = plot_paths(
+ "GBM",
+ qp.Lattice(2**8, seed=42),
+ t_final=5,
+ initial_value=50,
+ drift=0.1,
+ diffusion=0.2,
+ n=32,
+)
+```
+
+The generated sample paths are plotted below. The four panels show, respectively:
+
+- BM with lattice sampler ($T=1$, $S_0=1$, $\mu=0$, $\sigma^2=1$, 16 paths),
+- GBM with lattice sampler ($T=1$, $S_0=1$, $\mu=0$, $\sigma^2=1$, 16 paths),
+- GBM with IID sampler ($T=5$, $S_0=50$, $\mu=0.1$, $\sigma^2=0.2$, 32 paths),
+- GBM with lattice sampler ($T=5$, $S_0=50$, $\mu=0.1$, $\sigma^2=0.2$, 32 paths).
+
+
+
+
+
+
+
+
+Comparison of sample paths for BM and GBM using different samplers.
+
+
+## QuantLib vs QMCPy Comparison
+
+In this section, we compare QMCPy's `GeometricBrownianMotion` implementation with the industry-standard QuantLib library [6] to validate its accuracy and performance. The numerical results are summarized in the following table.
+
+Both libraries produce statistically equivalent GBM simulations that match theoretical values. QMCPy typically runs 1.5 to 2 times faster due to vectorized operations, lazy loading, and optimized memory management. More importantly, it demonstrates superior numerical accuracy (lower mean absolute errors) with Sobol, lattice, and Halton samplers, making it useful for research and high-performance applications. QuantLib remains the industry standard for production systems that require comprehensive support for financial modeling and risk management.
+
+| Method | Sampler | Mean | Std Dev | Mean Absolute Error | Std Dev Error | Mean Time (s) | Std Dev (s) | Speedup |
+| --- | --- | ---: | ---: | ---: | ---: | ---: | ---: | ---: |
+| QuantLib | Sobol | 100.88887227 | 0.00000000 | 4.23823737 | 21.23743882 | 0.00009037 | 0.00000000 | - |
+| QMCPy | Sobol | 103.66893699 | 1.17417404 | 1.45817264 | 20.06326479 | 0.00143112 | 0.00000000 | 0.08401794 |
+
+```python
+import QuantLib as ql
+import numpy as np
+
+
+def generate_quantlib_paths(
+ initial_value: float,
+ mu: float,
+ sigma: float,
+ maturity: float,
+ n_steps: int,
+ n_paths: int,
+ sampler_type: str = "IIDStdUniform",
+ seed: int = 7,
+) -> tuple:
+ """
+ Generate Geometric Brownian Motion paths using QuantLib.
+
+ Args:
+ initial_value: Initial value of the GBM process (S_0)
+ mu: Drift parameter
+ sigma: Volatility parameter (note: NOT diffusion coefficient)
+ maturity: Final time T
+ n_steps: Number of discretization time steps
+ n_paths: Number of paths to generate
+ sampler_type: Type of sampler ('IIDStdUniform' or 'Sobol')
+ seed: Random seed for reproducibility
+
+ Returns:
+ tuple: (paths, gbm) where paths has shape (n_paths, n_steps+1)
+ (includes initial value at t=0) and gbm is the
+ GeometricBrownianMotionProcess object
+
+ Raises:
+ ValueError: If sampler_type is not 'IIDStdUniform' or 'Sobol'
+ """
+ gbm = ql.GeometricBrownianMotionProcess(initial_value, mu, sigma)
+ times = ql.TimeGrid(maturity, n_steps)
+ dimension = n_steps
+ if sampler_type == "IIDStdUniform":
+ uniform_rng = ql.UniformRandomGenerator(seed)
+ sequence_gen = ql.GaussianRandomSequenceGenerator(
+ ql.UniformRandomSequenceGenerator(n_steps, uniform_rng)
+ )
+ path_gen = ql.GaussianPathGenerator(gbm, maturity, n_steps, sequence_gen, False)
+ paths = np.zeros((n_paths, n_steps + 1))
+ for i in range(n_paths):
+ sample_path = path_gen.next().value()
+ paths[i, :] = np.array([sample_path[j] for j in range(n_steps + 1)])
+ return paths, gbm
+ elif sampler_type == "Sobol":
+ uniform_rsg = ql.UniformLowDiscrepancySequenceGenerator(dimension, seed)
+ sequence_gen = ql.GaussianLowDiscrepancySequenceGenerator(uniform_rsg)
+ path_gen = ql.GaussianSobolMultiPathGenerator(
+ gbm, list(times), sequence_gen, False
+ )
+ paths = np.zeros((n_paths, n_steps + 1))
+ paths[:, 0] = initial_value # Set initial value
+ for i in range(n_paths):
+ sample_path = path_gen.next().value()
+ # For 1D process, get the first (and only) path
+ path_values = sample_path[0]
+ paths[i, :] = np.array([path_values[j] for j in range(n_steps + 1)])
+ return paths, gbm
+ else:
+ raise ValueError(
+ f"Unsupported sampler type: {sampler_type}. Use 'IIDStdUniform' or 'Sobol'"
+ )
+```
+
+In the next figure, the top row shows sample paths: QMCPy on the left
+and QuantLib on the right. The bottom-left panel overlays the marginal
+distribution at $t=1$, where both libraries yield nearly identical
+shapes with means and standard deviations matching the table above. The
+bottom-right panel is a QMCPy covariance heatmap with time 0 at the top
+of the $y$-axis; the variance increases along the diagonal with time and
+the off-diagonal structure follows $\min(t_i,t_j)$, consistent with the
+analytic form and the numerical matrix above.
+
+
+
+ QMCPy vs QuantLib comparison. Top: sample paths from QMCPy (left) and QuantLib (right). Bottom left: marginal distribution at $t=1$. Bottom right: QMCPy covariance heatmap.
+
+
+The evaluation of computational efficiency was done by creating comprehensive performance benchmarks comparing QMCPy and QuantLib across two key scaling dimensions. The benchmarks were performed using the `perfplot` library, which automatically handles warm-up, multiple runs, and statistical analysis to ensure reliable timing measurements.
+
+The following figure presents the results of our performance analysis. The left panel shows how execution time scales with the number of time steps while keeping the number of paths fixed at 4,096. Both libraries exhibit approximately linear scaling, but QMCPy demonstrates superior performance at smaller time step counts, with QuantLib becoming more competitive as the number of time steps increases. The right panel examines scaling behavior with respect to the number of paths while fixing the time steps at 252, representing a typical trading year. Here, QMCPy maintains a consistent performance advantage across all path counts, with the gap becoming more pronounced at higher path numbers. This performance difference is particularly relevant for Monte Carlo applications that require large numbers of simulation paths for accurate estimation.
+
+
+
+ GBM Performance Comparison: QuantLib vs QMCPy. Left plot shows performance scaling with number of time steps for fixed paths. Right plot shows performance scaling with number of paths for fixed number of time steps.
+
+
+To further validate these performance findings, we conducted an extended parameter sweep analysis across a broader range of configurations. The next figure presents the comprehensive results of this analysis, systematically examining performance across varying time steps and path counts.
+
+- Regarding accuracy, QMCPy's Sobol sampler generally achieves the lowest
+ mean absolute error (MAE), particularly when a larger number of paths
+ are used, reaching errors below $10^{-3}$. The lattice sampler also
+ provides competitive accuracy, especially with higher path counts.
+ While QuantLib's standard uniform IID sampler yields slightly lower
+ MAE for a larger number of paths compared to QMCPy's IID sampler, it
+ is notably slower.
+- Regarding speed, low-discrepancy samplers like Sobol, lattice, and
+Halton demonstrate superior convergence rates with increasing path counts compared to IID methods. With QMCPy, Sobol and lattice samplers offer the best speed-accuracy trade-off. QuantLib achieves comparable runtime to QMCPy's faster samplers, but without the accuracy benefits. The Halton sampler, while yielding the most accurate results, incurs significantly higher computational costs. These results highlight QMCPy's quasi-Monte Carlo methods as particularly well-suited for applications requiring high accuracy, with Sobol and lattice samplers providing an optimal balance of speed and precision for most practical scenarios.
+
+
+
+ Comprehensive parameter sweep performance analysis comparing QuantLib and QMCPy across varying time steps and path counts. The results demonstrate QMCPy's consistent performance advantages and superior scaling characteristics for both dimensions of the parameter space.
+
+
+## Internals
+
+The `GeometricBrownianMotion` class in QMCPy is engineered for speed, robustness, and mathematical correctness. Its design leverages object-oriented inheritance and vectorized operations, resulting in both flexibility and high performance. `GeometricBrownianMotion` inherits from `BrownianMotion`, which itself inherits from `Gaussian`. This layered design allows the GBM class to reuse and extend efficient implementations for Gaussian random vectors and BM increments. The constructor rigorously checks input parameters (e.g., positivity of initial value and diffusion, valid decomposition type), ensuring mathematical integrity and preventing run-time errors.
+
+The class uses vectorized NumPy operations to generate entire arrays of GBM paths in a single call, minimizing Python loops and maximizing computational throughput. Sample generation proceeds in two stages:
+
+1. The parent class `BrownianMotion` generates standard BM sample paths
+using the specified sampler (e.g., low-discrepancy lattice, IID uniform), with drift and diffusion handled in the mean and covariance structure.
+2. The GBM class transforms the BM samples via the exponential mapping
+above, performed in a fully vectorized fashion, ensuring that thousands of paths can be efficiently simulated.
+
+The class computes and stores the theoretical mean and covariance matrices for GBM at initialization, which can be used for validation and theoretical comparisons. Both mean and covariance are calculated using analytical formulas, leveraging [NumPy broadcasting](https://numpy.org/devdocs/user/basics.broadcasting.html) for efficient computations. Briefly, broadcasting in NumPy allows arithmetic operations between arrays of different shapes by automatically expanding the smaller array to match the shape of the larger array.
+
+The Gaussian and BM classes both implement Cholesky and PCA factorization of the covariance matrix
+
+$$
+\Sigma = L L^\top
+\quad\text{and}\quad
+\Sigma = P D P^\top,
+$$
+
+respectively. In the Cholesky method, one computes the lower-triangular
+$L=\operatorname{chol}(\Sigma)$ and obtains correlated increments via
+$X=LZ$, where $Z\sim\mathcal{N}(0,I)$. In the PCA approach, one first
+diagonalizes $\Sigma=PDP^\top$, then forms $X = P D^{1/2} Z$. By
+default PCA is used for its superior numerical stability in high
+dimensions and slightly lower cost when many eigenvalues are near zero.
+These correlated normals feed directly into the GBM update
+
+$$
+S_{t+\Delta t}
+= S_t \exp\!\Bigl((\mu - \tfrac12\sigma^2)\Delta t
++ \sigma\sqrt{\Delta t}\,X\Bigr),
+$$
+
+ensuring that the simulated paths respect the intended covariance structure and remain strictly positive, with strict-positivity checks raising warnings or errors if violated.
+
+## Conclusions and Future Work
+
+This blog demonstrates that QMCPy's quasi-Monte Carlo implementations provide significant advantages over traditional Monte Carlo methods for modeling geometric Brownian motion. QMCPy's approach combines superior numerical accuracy with enhanced computational efficiency, making it particularly well-suited for high-performance financial modeling applications.
+
+In the future, we aim to investigate ways to improve the runtime of the Halton sampler. For example, we may consider starting with Halton sampling points for high accuracy in early iterations, then switch to Sobol for faster convergence as sample size increases. It would also be interesting to experiment with ensemble sampling by running multiple samplers in parallel and combine the results using weighted averaging based on their relative accuracies.
+
+::: {.callout-tip}
+## Takeaways
+
+To the best of our knowledge, this blog presents the first publicly
+available benchmark comparing QuantLib and QMCPy for Geometric
+Brownian Motion simulations. Key takeaways include:
+
+- **Theoretical Validation**: QMCPy's GBM implementation correctly
+ preserves log-normality and matches theoretical moments, providing
+ reliable foundations for financial modeling applications.
+- **QMC Advantage**: Quasi-Monte Carlo methods (Sobol, lattice,
+ Halton) demonstrate superior convergence and accuracy compared to
+ traditional Monte Carlo, with Halton achieving the lowest mean
+ error and Sobol providing consistent performance across scenarios.
+- **Flexibility and Ease of Use**: Change one argument in QMCPy to
+ swap samplers or covariance decompositions, enabling rapid
+ experimentation and method comparison.
+- **Performance**: Vectorized operations and optimized memory
+ management produce 1.5-2x faster path generation than QuantLib.
+- **Robustness**: Superior numerical stability through PCA
+ decomposition and comprehensive error handling make QMCPy suitable
+ for both research and practical financial and scientific
+ applications.
+:::
+
+## References
+
+1. Glasserman, P. *Monte Carlo Methods in Financial Engineering*.
+ Springer-Verlag, New York, 2004.
+2. Hull, J. *Options, Futures, and Other Derivatives* (10th ed.).
+ Pearson, 2017.
+3. Ross, S. M. *Introduction to Probability Models* (11th ed.).
+ Academic Press, 2014.
+4. Choi, S.-C. T., Hickernell, F. J., Jagadeeswaran, R., McCourt, M. J.,
+& Sorokin, A. G. Quasi-Monte Carlo Software. In A. Keller (Ed.), *Monte Carlo and Quasi-Monte Carlo Methods*, 2022.
+5. Choi, S.-C. T., Hickernell, F., McCourt, M., & Sorokin, A. QMCPy:
+A quasi-Monte Carlo Python Library. [https://qmcsoftware.github.io/QMCSoftware/](https://qmcsoftware.github.io/QMCSoftware/). 2020.
+6. The QuantLib contributors. QuantLib: A free/open-source library for
+quantitative finance. Version 1.38. DOI: [10.5281/zenodo.1440997](https://doi.org/10.5281/zenodo.1440997). 2003--2025. [https://www.quantlib.org](https://www.quantlib.org).
+
+## Acknowledgments
+
+The authors thank Fred Hickernell, Joshua Jay Herman and Jiangrui Kang for their insightful feedback and help with the blog post.
diff --git a/blog/index.qmd b/blog/index.qmd
index abc05ad..bb72c75 100644
--- a/blog/index.qmd
+++ b/blog/index.qmd
@@ -18,6 +18,6 @@ listing:
The QMCSoftware blog shares accessible introductions, project updates, and
technical perspectives from across the community.
-The migrated articles establish a pattern for bringing the earlier QMCPy blog
-archive into the umbrella site. Additional posts will be reviewed and migrated
-incrementally.
+This archive brings together the QMCPy articles previously maintained with the
+project documentation. Technical documentation and executable notebooks remain
+with QMCPy and are linked from the relevant articles.
diff --git a/blog/linear-matrix-scrambling-and-digital-shift-for-halton/figures/ds.png b/blog/linear-matrix-scrambling-and-digital-shift-for-halton/figures/ds.png
new file mode 100644
index 0000000..ec4d264
Binary files /dev/null and b/blog/linear-matrix-scrambling-and-digital-shift-for-halton/figures/ds.png differ
diff --git a/blog/linear-matrix-scrambling-and-digital-shift-for-halton/figures/lms-ds.png b/blog/linear-matrix-scrambling-and-digital-shift-for-halton/figures/lms-ds.png
new file mode 100644
index 0000000..bff99b0
Binary files /dev/null and b/blog/linear-matrix-scrambling-and-digital-shift-for-halton/figures/lms-ds.png differ
diff --git a/blog/linear-matrix-scrambling-and-digital-shift-for-halton/figures/lms.png b/blog/linear-matrix-scrambling-and-digital-shift-for-halton/figures/lms.png
new file mode 100644
index 0000000..cfcda6c
Binary files /dev/null and b/blog/linear-matrix-scrambling-and-digital-shift-for-halton/figures/lms.png differ
diff --git a/blog/linear-matrix-scrambling-and-digital-shift-for-halton/index.qmd b/blog/linear-matrix-scrambling-and-digital-shift-for-halton/index.qmd
new file mode 100644
index 0000000..10de47f
--- /dev/null
+++ b/blog/linear-matrix-scrambling-and-digital-shift-for-halton/index.qmd
@@ -0,0 +1,152 @@
+---
+title: "Linear Matrix Scrambling and Digital Shift for Halton"
+author: "Aadit Jain"
+date: 2025-10-01
+date-format: "MMMM D, YYYY"
+description: "Linear matrix scrambling and digital shifts for randomized Halton sequences."
+categories:
+ - "QMCPy"
+ - "Halton Sequences"
+ - "Randomization"
+image: figures/lms-ds.png
+---
+
+This post introduces linear matrix scrambling and digital shifts for Halton sequences and shows how these randomizations affect point generation.
+
+## Introduction: Halton Sequences and Their Randomizations
+
+The Halton sequence is a common low-discrepancy sequence used for quasi-Monte Carlo simulations. It is based on the principle of using prime numbers as bases for each dimension. For example, base 2 is used for dimension 1, base 3 for dimension 2, base 5 for dimension 3, and so on. For each dimension, the index is converted from base 10 to the corresponding base, the digits are reversed, a decimal point is added before the reversed digits, and the result is converted back to base 10 to generate the sample coordinate.
+
+Here is an example illustrating how Halton samples are generated. For a 3-dimensional Halton object, the 29th index, which gives the 30th sample when counting from 1, is calculated as follows:
+
+$$
+i = 29 = 11101_{2} = 1002_{3} = 104_{5}.
+$$
+
+$$
+11101_{2} \to .10111_{2} = \frac{23}{32}.
+$$
+
+$$
+1002_{3} \to .2001_{3} = \frac{55}{81}.
+$$
+
+$$
+104_{5} \to .401_{5} = \frac{101}{125}.
+$$
+
+Hence, the 30th sample is
+$\left(\frac{23}{32}, \frac{55}{81}, \frac{101}{125}\right)$.
+
+Like lattice and Sobol sequences, Halton sequences have certain randomizations available to generate the samples. Previously, two randomizations of Halton implemented in QMCPy were `QRNG` [1] and `OWEN` [2]. `QRNG` uses optimized fixed permutations of the digits and also adds a random digital shift to the digits, while `OWEN` uses independent random permutations of the digits.
+
+This blog discusses the implementation of three new randomizations for Halton: linear matrix scrambling (LMS), digital shift (DS), and linear matrix scrambling plus digital shift (LMS_DS). These randomizations are commonly used for digital nets but have only recently been explored for Halton sequences. They help provide unbiased estimates for QMC integration, and variance under a linear matrix scramble can be better than variance under only a random digital shift.
+
+## LMS, DS, and LMS_DS
+
+1. `LMS`: Linear matrix scrambling of Halton [3]. Based on the bases, a
+different scrambling matrix is generated for each dimension. The lower triangle is random between 0 and base minus 1, the diagonal is random between 1 and base minus 1, and the upper triangle is all zeros. After the indexes are converted to their base representations and a decimal point has been added before the reversed digits, their dot product is computed with the scrambling matrix. After computing the dot product, the scrambled indexes or coefficients are converted to base 10 to generate the samples.
+2. `DS`: Digital shift. Based on the bases, a different vector is
+generated for each dimension, with entries random between 0 and base minus 1. After the indexes are converted to their base representations and a decimal point has been added before the reversed digits, they are added to the vector and then converted to base 10 to generate the samples.
+3. `LMS_DS`: A combination of linear matrix scrambling and digital shift.
+Linear matrix scrambling of Halton is applied first; then, before converting to base 10, the digital shift is applied to the scrambled indexes or coefficients.
+
+`LMS` includes the origin as the first point in the sequence, but `DS` prevents this. This makes `LMS_DS` the recommended randomized option among these three.
+
+## Plot Examples of LMS, DS, and LMS_DS
+
+```python
+import qmcpy as qp
+
+dimension = 3
+
+lms_halton = qp.Halton(dimension, randomize="LMS")
+ds_halton = qp.Halton(dimension, randomize="DS")
+lms_ds_halton = qp.Halton(dimension, randomize="LMS_DS")
+
+fig1, ax1 = qp.plot_proj(
+ lms_halton,
+ n=2**5,
+ d_horizontal=range(dimension),
+ d_vertical=range(dimension),
+ math_ind=False,
+)
+fig2, ax2 = qp.plot_proj(
+ ds_halton,
+ n=2**5,
+ d_horizontal=range(dimension),
+ d_vertical=range(dimension),
+ math_ind=False,
+)
+fig3, ax3 = qp.plot_proj(
+ lms_ds_halton,
+ n=2**5,
+ d_horizontal=range(dimension),
+ d_vertical=range(dimension),
+ math_ind=False,
+)
+
+fig1.suptitle("LMS")
+fig2.suptitle("DS")
+fig3.suptitle("LMS_DS")
+```
+
+
+
+ Figure 1: LMS Halton plot.
+
+
+
+
+ Figure 2: DS Halton plot.
+
+
+
+
+ Figure 3: LMS_DS Halton plot.
+
+
+More plot examples can be seen in the [Linear Matrix Scrambling and Digital Shift for Halton notebook](https://github.com/QMCSoftware/QMCSoftware/blob/develop/demos/linear-scrambled-halton.ipynb).
+
+## Speed Comparison Between Halton Randomization Methods
+
+```python
+import qmcpy as qp
+from time import time
+
+dimension = 25
+samples = 100000
+rand_options = ["QRNG", "OWEN", "LMS", "DS", "LMS_DS"]
+
+for rand_option in rand_options:
+ t_start = time()
+ qp.Halton(dimension, randomize=rand_option).gen_samples(samples, warn=False)
+ t_end = time()
+ print(f"Time to generate samples for {rand_option}= {t_end - t_start}")
+```
+
+Output from the original timing run:
+
+```text
+Time to generate samples for QRNG= 0.355226993560791
+Time to generate samples for OWEN= 2.15480637550354
+Time to generate samples for LMS= 1.348050832748413
+Time to generate samples for DS= 1.0668678283691406
+Time to generate samples for LMS_DS= 2.5354738235473633
+```
+
+Through this speed comparison, we can see that `LMS_DS` is slower than the other randomization techniques.
+
+## Conclusion
+
+`LMS`, `DS`, and `LMS_DS` provide newer ways to randomize Halton points and obtain a wider variety of samples. The digital shift prevents unrandomized Halton and scrambled Halton from including the origin as the first point, which allows these sequences to be used with different true measure objects and integral approximations where the origin may be problematic.
+
+## References
+
+1. Hofert, M. & Lemieux, C. `qrng`: (Randomized) Quasi-Random Number
+Generators. R package version 0.0-7. 2019. [https://CRAN.R-project.org/package=qrng](https://CRAN.R-project.org/package=qrng).
+2. Owen, A. B. A randomized Halton algorithm in R. 2017.
+ [arXiv:1706.02808](https://arxiv.org/abs/1706.02808) [stat.CO].
+3. Owen, A. B. & Pan, Z. Gain coefficients for scrambled Halton points.
+ 2023. [arXiv:2308.08035](https://arxiv.org/abs/2308.08035)
+ [math.NA].
diff --git a/blog/qei-with-qmcpy/figures/qei_convergence.png b/blog/qei-with-qmcpy/figures/qei_convergence.png
new file mode 100644
index 0000000..121d6a2
Binary files /dev/null and b/blog/qei-with-qmcpy/figures/qei_convergence.png differ
diff --git a/blog/qei-with-qmcpy/figures/test_problem.png b/blog/qei-with-qmcpy/figures/test_problem.png
new file mode 100644
index 0000000..2862498
Binary files /dev/null and b/blog/qei-with-qmcpy/figures/test_problem.png differ
diff --git a/blog/qei-with-qmcpy/index.qmd b/blog/qei-with-qmcpy/index.qmd
new file mode 100644
index 0000000..ec22fe3
--- /dev/null
+++ b/blog/qei-with-qmcpy/index.qmd
@@ -0,0 +1,84 @@
+---
+title: "qEI with QMCPy"
+author: "Michael McCourt"
+date: 2020-07-27
+date-format: "MMMM D, YYYY"
+description: "Using QMCPy low-discrepancy samples to estimate q-Expected Improvement in Bayesian optimization."
+categories:
+ - "QMCPy"
+ - "Bayesian Optimization"
+ - "Quasi-Monte Carlo"
+image: figures/test_problem.png
+---
+
+This post demonstrates how QMCPy low-discrepancy samples can improve Monte Carlo estimation of q-Expected Improvement in Bayesian optimization.
+
+Quasi-Monte Carlo methods (QMC) are a valuable tool for sampling random variables in a structured fashion. This allows for computing key statistics of random variables [more efficiently](../why-add-q-to-mc/index.qmd) than with [i.i.d. sampling](https://en.wikipedia.org/wiki/Independent_and_identically_distributed_random_variables). Such quantities can play fundamental roles in larger algorithms, making their efficient computation fundamental to practical implementations of numerous applications. This was the motivation for creating the QMCPy library. In this post, we demonstrate the use of QMC methods in computing a key quantity in Bayesian optimization.
+
+## Bayesian Optimization
+
+[Bayesian optimization](https://arxiv.org/abs/1807.02811), also called many other names, including [sequential model based optimization](https://link.springer.com/chapter/10.1007/978-3-642-25566-3_40) or [Gaussian process optimization](https://icml.cc/Conferences/2010/papers/422.pdf), is a broad class of algorithms which involves alternately building statistical models of scattered data and optimally sampling for further data to try to optimize a given function. To conduct this optimal sampling process, a strategy must be applied to take the statistical model and appropriately balance a desire to explore the optimization domain against a desire to exploit the high performing values found thus far and search in their proximity. We refer to this strategy as the [acquisition function](https://www.cse.wustl.edu/~garnett/cse515t/spring_2015/files/lecture_notes/12.pdf).
+
+Probably the most popular acquisition function is
+[Expected Improvement](https://link.springer.com/article/10.1023/A:1008306431147)
+(EI), which is widely used both because of its closed form for a
+[Gaussian process](https://sigopt.com/blog/intuition-behind-gaussian-processes/)
+surrogate model and because of its effectiveness. When trying to push
+beyond the sequential model of Bayesian optimization to allow for
+parallelism, the EI acquisition function can be redefined for optimally
+choosing the next $q$ points at which to sample. This qEI acquisition
+function, first defined
+[here](https://link.springer.com/chapter/10.1007/978-3-642-10701-6_6),
+with more recent discussion
+[here](https://pubsonline.informs.org/doi/abs/10.1287/opre.2019.1966),
+unfortunately does not have a simple closed form, and is generally
+computed through Monte Carlo estimation.
+
+## qEI Estimation
+
+We consider a simple problem to demonstrate the value of QMC in estimating this qEI quantity. The left panel of Figure 1 depicts 5 noisy observations of a 1D function. We fit a Gaussian process (GP) model to this data, and provide examples of posterior draws in the center panel.
+
+
+
+ Figure 1: left: 5 observed points, drawn from a function, with error bars displayed. center: 128 quasi-random posterior draws from a Gaussian process fit to the observed points. right: The qEI associated with $q=2$ points given this Gaussian process and the observed points; note the symmetry across the line $y=x$, which occurs because the points will be sampled at the same time.
+
+
+The right panel of Figure 1 depicts the qEI quantity for $q=2$ future
+points to be sampled. qEI is defined, for a maximization problem, with a
+$q$-dimensional integral as
+
+$$
+\operatorname{qEI}(x_1, \ldots, x_q)
+=
+\int_{\mathbb{R}^q}
+\max_{1 \le i \le q}(y_i - y^*)_+
+p_{Y_{x_1, \ldots, x_q}}(y_1, \ldots, y_q)
+\, \mathrm{d}y_1 \cdots \mathrm{d}y_q,
+$$
+
+where $Y_{x_1, \ldots, x_q}$ is the joint posterior GP distribution at
+$x_1, \ldots, x_q$ and $y^*$ is the best value observed thus far. Note
+the symmetric nature of the qEI: since both of the $q=2$ points will be
+evaluated simultaneously, the designation of "first" and "second" is
+arbitrary. The highest qEI value occurs when sampling near $x_1=0$ and
+$x_2=0.6$, points which are near the highest valued draws of the
+posterior.
+
+This integral is generally estimated using Monte Carlo, as discussed
+[here](http://www.cs.ubc.ca/labs/beta/EARG/stack/2010_CI_Ginsbourger-ParallelKriging.pdf)
+and
+[here](http://papers.nips.cc/paper/8194-maximizing-acquisition-functions-for-bayesian-optimization),
+among other places, with i.i.d. draws. Using Quasi-Monte Carlo, we can
+reach the same level of accuracy with fewer GP posterior draws. Figure 2
+shows the results of using different quasi-random sequences available in
+QMCPy to estimate qEI for $q=5$ next points, located at
+$0.158$, $0.416$, $0.465$, $0.718$, and $0.935$.
+
+
+
+ Figure 2: Monte Carlo estimation, computed both using QMCPy's i.i.d. sampler and NumPy's normal sampler, is compared to QMC methods available in QMCPy. The Sobol' and Lattice methods perform comparably, as expected, and converge at roughly $\mathcal{O}(N^{-1})$, in contrast to the $\mathcal{O}(N^{-1/2})$ convergence of i.i.d. sampling. Integrals were estimated 50 times with different random seeds: the medians (solid lines) and interquartile ranges (shaded regions) are plotted.
+
+
+The integrals present in Bayesian optimization are another example of a computation which benefits from QMC. More efficient acquisition function computation gives us the ability to conduct [more complicated strategies](http://www.auai.org/uai2020/proceedings/124_main_paper.pdf) in a feasible amount of time. Check out the [QMCPy documentation](https://qmcsoftware.github.io/QMCSoftware/) or [contact us](https://github.com/QMCSoftware/QMCSoftware/issues) to learn more about the QMCPy project and how QMC can help your work.
+
+For a related executable notebook, see the [qEI demo for blog](https://qmcsoftware.github.io/QMCSoftware/demos/qei-demo-for-blog/).
diff --git a/blog/random-lattice-generators-are-not-bad/figures/lattice-points.png b/blog/random-lattice-generators-are-not-bad/figures/lattice-points.png
new file mode 100644
index 0000000..cd2bfed
Binary files /dev/null and b/blog/random-lattice-generators-are-not-bad/figures/lattice-points.png differ
diff --git a/blog/random-lattice-generators-are-not-bad/figures/mean-vs-median.png b/blog/random-lattice-generators-are-not-bad/figures/mean-vs-median.png
new file mode 100644
index 0000000..cb065fc
Binary files /dev/null and b/blog/random-lattice-generators-are-not-bad/figures/mean-vs-median.png differ
diff --git a/blog/random-lattice-generators-are-not-bad/index.qmd b/blog/random-lattice-generators-are-not-bad/index.qmd
new file mode 100644
index 0000000..3031640
--- /dev/null
+++ b/blog/random-lattice-generators-are-not-bad/index.qmd
@@ -0,0 +1,217 @@
+---
+title: "Random Lattice Generators are Not Bad"
+author: "Bocheng David Zhang"
+date: 2024-05-13
+date-format: "MMMM D, YYYY"
+description: "Random lattice generating vectors and comparisons of mean, median, and randomly shifted lattice rules."
+categories:
+ - "QMCPy"
+ - "Lattice Rules"
+ - "Quasi-Monte Carlo"
+image: figures/lattice-points.png
+---
+
+This post discusses random lattice generating vectors in QMCPy and compares mean, median, and randomly shifted lattice rules.
+
+Generating vectors are used by lattice generators to compute point sets. Previous works [1, 2, 3, 4] commonly applied greedy component-by-component (CBC) algorithms to construct generating vectors. However, this process is dependent on weight vectors and the decay of Fourier coefficients. To this end, Takashi Goda and Pierre L'Ecuyer suggested in their work [5] that generating vectors do not require a predetermined weight vector and a measure of the decay of Fourier coefficients to reach high precision; instead, random generators can also achieve desirable results. Recently, we implemented random generating vectors into QMCPy, and the results were quite promising. This blog explores the usage of random generating vectors in QMCPy through code examples. Before that, however, we shall consider some mathematics behind generating vectors.
+
+## Mathematics of Generating Vectors
+
+Given $N > 2$ and generating vector
+$z \in \{1,\dots,N-1\}^d$, define
+$P_N := P_{N,z} := \{\boldsymbol{x}_0,\boldsymbol{x}_1,\dots,
+\boldsymbol{x}_{N-1}\} \subset [0,1)^d$, where
+
+$$
+x_n = \left\{\frac{nz}{N}\right\},
+\qquad n = 0,1,\dots,N-1,
+$$
+
+and $\{s\} = s-\lfloor s \rfloor$ denotes the fractional part of each
+component. The significance of the generating vector is that it
+determines the point set $P_N$.
+
+The following are plots of lattices generated by QMCPy when $N$ is a
+power of 2:
+
+
+
+ Figure 1: Lattices generated using generating_vector=16 and seed=136. These plots can be conditionally reproduced in the lattice random generator notebook.
+
+
+While samples in [5] were produced by a prime sample size, QMCPy
+utilizes sample sizes that are powers of 2 due to its expectation of
+extensible generating vectors. The usage of the generating vector comes
+from quasi-Monte Carlo integration using a point set $P_N$:
+
+$$
+I(f) := \int_{[0,1)^d} f(x) \, \mathrm{d}x
+\approx I(f;P_N) = \frac{1}{N}\sum_{n=0}^{N-1} f(x_n).
+$$
+
+In their work, Goda and L'Ecuyer approximated $I(f)$ by the median
+rule,
+$\operatorname{median}(I(f;P_{N,z_1}),\dots,I(f;P_{N,z_r}))$ [5],
+while QMCPy utilizes mean rules,
+$\operatorname{mean}(I(f;P_{N,z_1}),\dots,I(f;P_{N,z_r}))$. We
+performed numerical experiments in the final section to have a glimpse
+of the efficiency of both rules.
+
+## Code Examples
+
+In this section, we will explore the basic features of the `Lattice` class and the `gen_samples` method. For further documentation, see the [QMCPy Lattice documentation](https://qmcpy.readthedocs.io/en/latest/algorithms.html#module-qmcpy.discrete_distribution.lattice.lattice).
+
+The generating vector is the core of the `Lattice` object. Currently, QMCPy enables the following types of cubature schemes:
+
+1. A hard-coded $d$-dimensional array.
+2. A file that contains a hard-coded generating vector.
+3. A totally random generator produced by integer input.
+
+We will focus on the recently developed third type of generating vector because it is a direct application of the mathematics discussed above.
+
+## Lattice Declaration and the `gen_samples` Function
+
+A `Lattice` object in QMCPy requires the dimension and the generating vector of choice. Other arguments such as `randomize` or `seed` are optional.
+
+The following code is a short example used to illustrate the declaration of a `Lattice` object and the `gen_samples` function.
+
+```python
+import qmcpy as qp
+lattice = qp.Lattice(dimension=2, generating_vector=21, seed=120) # intialize the lattice
+print(lattice) # print information about the lattice
+print(lattice.gen_samples(n=4)) # print the first 4 points in the lattice
+```
+
+The output is listed below:
+
+```text
+Lattice (AbstractLDDiscreteDistribution)
+ d 2^(1)
+ replications 1
+ randomize SHIFT
+ gen_vec_source random
+ order RADICAL INVERSE
+ n_limit 2^(21)
+ entropy 120
+[[0.34548142 0.46736834]
+ [0.84548142 0.96736834]
+ [0.59548142 0.21736834]
+ [0.09548142 0.71736834]]
+```
+
+## Integration
+
+To integrate in QMCPy, one needs to declare the dimension, the
+tolerance, a low discrepancy sequence, and the true measure used to
+transform the objective function to the $d$-dimensional unit cube. In
+the following example, the Gaussian measure is applied using a lattice
+as the low discrepancy sequence. The Keister function [6] is used as an
+example.
+
+```python
+import qmcpy as qp
+
+d = 5
+tol = 1E-3
+
+data_random = qp.CubQMCLatticeG(
+ qp.Keister(
+ qp.Gaussian(
+ qp.Lattice(d, generating_vector=26),
+ mean=0, covariance=1/2)),
+ abs_tol = tol).integrate()[1]
+print("Integration data from a random lattice generator:")
+print(data_random)
+
+data_default = qp.CubQMCLatticeG(
+ qp.Keister(
+ qp.Gaussian(
+ qp.Lattice(d),
+ mean=0, covariance=1/2)),
+ abs_tol = tol).integrate()[1]
+print("\nIntegration data from the default lattice generator:")
+print(data_default)
+```
+
+Abbreviated output is listed below:
+
+```text
+Integration data from a random lattice generator:
+LDTransformData (AccumulateData Object)
+ solution 1.136
+ comb_bound_low 1.135
+ comb_bound_high 1.137
+ comb_flags 1
+ n_total 2^(17)
+ n 2^(17)
+ time_integrate 0.642
+Lattice (DiscreteDistribution Object)
+ d 5
+ dvec [0 1 2 3 4]
+ randomize 1
+ order natural
+ gen_vec [ 1 29661439 12472787 51447409 58451577]
+ entropy 39936265299936103070191134814877412899
+ spawn_key ()
+
+Integration data from the default lattice generator:
+LDTransformData (AccumulateData Object)
+ solution 1.135
+ comb_bound_low 1.134
+ comb_bound_high 1.135
+ comb_flags 1
+ n_total 2^(17)
+ n 2^(17)
+ time_integrate 0.626
+Lattice (DiscreteDistribution Object)
+ d 5
+ dvec [0 1 2 3 4]
+ randomize 1
+ order natural
+ gen_vec [ 1 182667 469891 498753 110745]
+ entropy 230872234427376376523997454058047592711
+ spawn_key ()
+```
+
+One can see that the default generator performs slightly better than the
+random generator, as the time to integrate is about $0.16$ seconds
+lower.
+
+## Comparison Between Lattice Generators
+
+To have a further glimpse into the performance of lattice generators, we
+compared the error of three types of random generators with respect to
+the sample size. The blue line depicts a random generator that applies
+the median rule
+$\operatorname{median}(I(f;P_{N,z_1}),\dots,I(f;P_{N,z_r}))$; the
+orange line depicts a generator that applies the mean rule
+$\operatorname{mean}(I(f;P_{N,z_1}),\dots,I(f;P_{N,z_r}))$; the green
+line is a randomly shifted hard-coded generator.
+
+
+
+ Figure 2: A comparison between lattice generators. This plot can be conditionally reproduced in the lattice random generator notebook.
+
+
+Here, we used sample sizes $N$ ranging from $2^6$ to $2^{18}$ and
+$r = 11$. We compared the results of integrating the 2-dimensional
+Keister integral over each sample size using each type of lattice
+generator. To reduce sampling variance, we repeated the trials $25$
+times and computed the averaged result.
+
+As shown in the plot, the mean of random shifts (green) outperforms the random generator using median rules (blue), which in turn outperforms the random generator using mean rules (orange). These results support findings in [5]. More numerical experiments under different circumstances should be conducted before making a conclusion, but current work suggests that random lattice generators have a lot of potential.
+
+## References
+
+1. Korobov, N. M. The approximate computation of multiple integrals.
+ *Doklady Akademii Nauk SSSR* 124, 1207-1210 (1959).
+2. Sloan, I. H. QMC integration: beating intractability by weighting
+the coordinate directions. In *Monte Carlo and Quasi-Monte Carlo Methods 2000* (eds. Fang, K. T., Hickernell, F. J., & Niederreiter, H.) 103-123 (Springer-Verlag, Berlin, 2002).
+3. Kuo, F. Y. Component-by-component constructions achieve the optimal
+rate of convergence for multivariate integration in weighted Korobov and Sobolev spaces. *Journal of Complexity* 19, 301-320 (2003).
+4. Nuyens, D. & Cools, R. Fast component-by-component construction. In
+*Monte Carlo and Quasi-Monte Carlo Methods 2004* (eds. Niederreiter, H. & Talay, D.) 373-387 (Springer-Verlag, Berlin, 2006).
+5. Goda, T. & L'Ecuyer, P. Construction-free median quasi-Monte Carlo
+rules for function spaces with unspecified smoothness and general weights. *SIAM Journal on Scientific Computing* 44, A2765-A2788 (2022). [https://doi.org/10.1137/22M1473625](https://doi.org/10.1137/22M1473625)
+6. Keister, B. D. Multidimensional quadrature algorithms. *Computers in
+ Physics* 10, 119-122 (1996).
diff --git a/blog/safe-handling-of-qmc-points/index.qmd b/blog/safe-handling-of-qmc-points/index.qmd
new file mode 100644
index 0000000..47aa023
--- /dev/null
+++ b/blog/safe-handling-of-qmc-points/index.qmd
@@ -0,0 +1,180 @@
+---
+title: "Safe Handling of QMC Points"
+author: "Art B. Owen"
+date: 2020-08-31
+date-format: "MMMM D, YYYY"
+description: "Why QMC points require more care than IID samples and why skipping, thinning, and arbitrary sample sizes can be unsafe."
+categories:
+ - "Quasi-Monte Carlo"
+ - "Low-Discrepancy Sequences"
+ - "Best Practices"
+---
+
+This post explains why QMC points should not be handled like IID samples and warns against unsafe practices such as skipping, thinning, and arbitrary sample sizes.
+
+A Quasi-Monte Carlo construction of $n$ points in $d$ dimensions may
+look like IID points, but they must be used with a bit more care.
+Because QMC can give errors that are $o(1/n)$ as $n \to \infty$,
+changing or ignoring even one point can change the estimate by an amount
+much larger than the error would have been and worsen the convergence
+rate. As a result, certain practices that fit quite naturally and
+intuitively with MC points are very detrimental to QMC performance.
+Operations like burn-in, thinning, and even using a round number sample
+size, like a power of ten, can degrade QMC effectiveness or even make it
+converge to the wrong answer. The safe way to use QMC points is to take
+all $n$ points produced, after applying a randomization to avoid
+singularities and to support uncertainty quantification.
+
+## Introduction
+
+This note arose from a discussion of quasi-Monte Carlo (QMC) and randomized quasi-Monte Carlo (RQMC) software during and following the plenary tutorial at MCQMC 2020 by Fred Hickernell. Common ways of handling IID points can fail to work for (R)QMC points. A longer discussion of this point is available at [arXiv:2008.08051](https://arxiv.org/abs/2008.08051).
+
+QMC sampling methods provide a set of $n$ points in $[0,1]^d$ that we
+can use instead of a sample of $\mathcal{U}[0,1]^d$ points. We can apply
+transformations to them to simulate non-uniform distributions and
+domains other than the unit cube. Then the resulting points can be used
+to estimate an expectation or just to explore the input to a function.
+
+If the points are
+$\boldsymbol{x}_1, \dots, \boldsymbol{x}_n \in [0,1]^d$, we may
+estimate
+$\mu = \int_{[0,1]^d} f(\boldsymbol{x}) \, \mathrm{d}\boldsymbol{x}$
+by
+
+$$
+\hat{\mu}
+= \frac{1}{n} \sum_{i=1}^n f(\boldsymbol{x}_i),
+$$
+
+just as we would have done with
+$\boldsymbol{x}_i \overset{\text{iid}}{\sim} \mathcal{U}[0,1]^d$. The
+function $f(\cdot)$ subsumes transformations as well as the integrand of
+interest in the transformed space.
+
+Plain QMC points are deterministic. Randomizing them in one of several
+possible ways makes them individually uniformly distributed while
+preserving the low discrepancy structure that makes them valuable for
+integration. The resulting RQMC methods allow uncertainty quantification
+via replication. If it is important to be accurate, then it must also be
+important to know that you were accurate and to show that you were
+accurate. A plain $t$-test based confidence interval, or better yet a
+bootstrap-$t$ confidence interval for $\mu$, then lets one estimate
+accuracy.
+
+Bootstrap-$t$ works very well even with a modest number of replicates.
+We might want a modest number $R$ of replicates because the root mean
+squared error (RMSE) decreases proportionally to $1/\sqrt{R}$ as the
+number of replicates increases but often faster than $1/\sqrt{n}$ as
+the number of sample points increases. The work involved is proportional
+to $nR$.
+
+A second reason to randomize is that QMC points are really designed for
+Riemann integrable functions. Those are necessarily bounded. If
+$\hat{\mu} \to \mu$ whenever the star discrepancy of
+$\boldsymbol{x}_1, \dots, \boldsymbol{x}_n$ converges to zero, then it
+must hold that $f$ is Riemann integrable. That is, if $f$ is not Riemann
+integrable, as for instance it would be if it were unbounded, then there
+are sequences of inputs with vanishing star discrepancy for which
+$\hat{\mu} - \mu$ does not converge to zero.
+
+It is safer to randomize. Nested uniform scrambles ensure that
+$\hat{\mu} \to \mu$ with probability one under the weak condition that
+$f \in L^{1+\epsilon}[0,1]^d$ for some $\epsilon > 0$. That is,
+
+$$
+\int_{[0,1]^d}
+\left|f(\boldsymbol{x})\right|^{1+\epsilon}
+\, \mathrm{d}\boldsymbol{x}
+< \infty,
+$$
+
+and $f$ is measurable.
+
+Because (R)QMC points look so similar to plain IID points, many users and software implementations handle (R)QMC points in inefficient or even unsafe ways that would be no problem for IID points.
+
+## Sample Sizes
+
+(R)QMC points are usually constructed as a finite sequence of points for
+a specific sample size $n$, such as $n=2^m$ or $n=p$ for a large prime
+number $p$. If one uses only a round number such as $1000$ of them, then
+the result will ordinarily be much less effective than using them all
+and can possibly even fail to sample a portion of $[0,1]^d$. Those
+$1000$ points might easily be less effective than using a smaller
+sequence of $512$ points. As for antibiotics, one should use the whole
+sequence.
+
+## Skipping or Burn-In
+
+For IID points, we do as well with
+$\boldsymbol{x}_{B+1}, \dots, \boldsymbol{x}_{B+n}$ for any $B \ge 0$.
+Taking $B > 0$ is a kind of burn-in that actually has an advantage in
+Markov chain Monte Carlo, where the points may only approach their
+desired distribution. For RQMC points, skipping even one observation can
+make the rate of convergence worse. In the case of scrambled nets,
+taking $B=1$ can turn the RMSE from approximately
+$\mathcal{O}(n^{-3/2})$ to approximately $\mathcal{O}(n^{-1})$.
+
+The reason that people often skip the first point is that this first
+point is often equal to $(0,0,\dots,0)$. Such a point is then
+problematic when $f$ maps $[0,1]^d$ onto $\mathbb{R}^d$, as it would
+when using a transformation to induce a Gaussian distribution before
+evaluating the quantity of interest. The point at the origin can map to
+an infinite point or even result in `not a number`.
+
+If one uses RQMC, then that first point ends up with the
+$\mathcal{U}[0,1]^d$ distribution, as do all the others. That avoids the
+problem of singularities at least mathematically. One might still hit a
+singularity in a floating point representation if one is extremely
+unfortunate. That possibility is also there with QMC, and plain QMC does
+not have the same assurance of avoiding singularities that RQMC has.
+
+## Thinning
+
+In MCMC one often takes every $k$-th point for reasons of storage or
+computational efficiency. In IID sampling taking every $k$-th point
+would be statistically equivalent to taking an equal number of
+consecutive points. If we use $\boldsymbol{x}_{ki}$ for integer $k > 1$
+and $i=1, \dots, n$ in (R)QMC, the result can be disastrously bad. For
+instance, the van der Corput sequence in $[0,1]$ alternates between
+values in $[0,1/2)$ and values in $[1/2,1)$. Taking every second point
+would ignore half of the domain. The first component of a Sobol'
+sequence is ordinarily the van der Corput sequence.
+
+Thinning (R)QMC points can be extremely dangerous. It should not be done without some very careful mathematical explanation of why it might be ok in some special setting.
+
+## van der Corput Sequences
+
+These are for $d=1$, so $x_i \in [0,1]$. Any consecutive $2^m$ points
+of the van der Corput sequence are a digital net and hence have some
+good discrepancy properties. The same holds for generalizations of van
+der Corput to bases $b > 2$. There any $b^m$ consecutive points are a
+digital net. So van der Corput points are an exception. If we use
+burn-in, we still get a digital net and thus still get low discrepancy.
+We should take care of the chosen sample size, preferring $n$ to be a
+power of $b$. If the powers of $b$ are too far apart for our purposes,
+then an integer multiple of a power of $b$ is next best. That only makes
+a difference when $b > 2$.
+
+We should not thin van der Corput sequences.
+
+## Halton Sequences
+
+Halton sequences are somewhat robust to burn-in and using round numbers.
+Each of the $d$ component variables of a Halton sequence is a van der
+Corput sequence in a different base. Usually the first $d$ prime numbers
+are used.
+
+For modestly large $d$, the special values of $n$ are so large and so
+far apart that we can consider that there simply are no specially good
+sample sizes. Think of making $n$ divisible by a power of the product of
+the first $d$ prime numbers. Even the first such value may be too large
+to use. When no feasible sample sizes are very good, then maybe there is
+no particular harm from using a power of ten.
+
+Halton sequences start at the origin, which is problematic as described
+above. We can easily skip that point in Halton sequences because there
+are no especially good ranges. It may even be advantageous to use a very
+large burn-in for the Halton sequence because the initial points for
+large $d$ have unpleasant striping artifacts.
+
+It is however safer to randomize the Halton sequence. Scrambling the Halton sequence counters those striping artifacts more surely than a burn-in would. It also moves the point at the origin to a uniformly distributed random point. This is another instance where RQMC is safer and more effective than plain QMC.
diff --git a/blog/scipywrapper/figures/fig01_mvn_indep_vs_dep.png b/blog/scipywrapper/figures/fig01_mvn_indep_vs_dep.png
new file mode 100644
index 0000000..ac87845
Binary files /dev/null and b/blog/scipywrapper/figures/fig01_mvn_indep_vs_dep.png differ
diff --git a/blog/scipywrapper/figures/fig02_zero_inflated_joint.png b/blog/scipywrapper/figures/fig02_zero_inflated_joint.png
new file mode 100644
index 0000000..100d9e2
Binary files /dev/null and b/blog/scipywrapper/figures/fig02_zero_inflated_joint.png differ
diff --git a/blog/scipywrapper/figures/fig03_accept_reject_mc_vs_qmc.png b/blog/scipywrapper/figures/fig03_accept_reject_mc_vs_qmc.png
new file mode 100644
index 0000000..c65e43b
Binary files /dev/null and b/blog/scipywrapper/figures/fig03_accept_reject_mc_vs_qmc.png differ
diff --git a/blog/scipywrapper/figures/fig04_custom_triangular_marginal.png b/blog/scipywrapper/figures/fig04_custom_triangular_marginal.png
new file mode 100644
index 0000000..f2e70a8
Binary files /dev/null and b/blog/scipywrapper/figures/fig04_custom_triangular_marginal.png differ
diff --git a/blog/scipywrapper/figures/fig05_multivariate_student_t_joint.png b/blog/scipywrapper/figures/fig05_multivariate_student_t_joint.png
new file mode 100644
index 0000000..38f1703
Binary files /dev/null and b/blog/scipywrapper/figures/fig05_multivariate_student_t_joint.png differ
diff --git a/blog/scipywrapper/index.qmd b/blog/scipywrapper/index.qmd
new file mode 100644
index 0000000..5199206
--- /dev/null
+++ b/blog/scipywrapper/index.qmd
@@ -0,0 +1,181 @@
+---
+title: "Extending QMCPy's `SciPyWrapper` to Support Dependent and Custom Distributions"
+author: "Laasya Priya Vemuri"
+date: 2026-04-18
+date-format: "MMMM D, YYYY"
+description: "Extending SciPyWrapper for dependent joint distributions, custom marginals, and validation warnings."
+categories:
+ - "QMCPy"
+ - "Probability Distributions"
+ - "SciPy"
+image: figures/fig01_mvn_indep_vs_dep.png
+---
+
+This post describes the updated `SciPyWrapper` interface for dependent joint distributions, custom marginals, and validation warnings.
+
+#### Introduction
+
+`SciPyWrapper` is one of QMCPy's most useful interfaces: it maps low-discrepancy samples on `[0,1]^d` to probability models used in simulation and integration. The original interface served independent univariate marginals from `scipy.stats` effectively, but many modern workflows require richer structure, including dependence and custom distribution families.
+
+This work extends `qmcpy.true_measure.SciPyWrapper` in three directions:
+
+1. Support for dependent joint distributions.
+2. Support for custom univariate distributions with SciPy-like behavior.
+3. Built-in warnings for clearly problematic custom distribution definitions.
+
+#### Motivation
+
+Applied uncertainty quantification rarely remains fully independent. In practice, modelers often need to represent correlated random variables, mixture behavior with atoms (such as zero inflation), and distribution families that are not directly available in SciPy. Prior to this extension, handling those cases required workarounds outside the normal `SciPyWrapper` flow.
+
+The goal of this work was to remove that gap while preserving backward compatibility for existing code.
+
+#### Methodology
+
+The updated implementation keeps the original independent-marginal path unchanged and adds a second path for dependent joint models.
+
+##### 1. Joint-distribution mode
+
+In joint mode, `SciPyWrapper` accepts an object with:
+
+1. `transform(u)` that maps `u in [0,1]^d` into samples in `R^d`.
+2. A dimension descriptor (`dim` or `dimension`) for consistency checks.
+3. Optional `logpdf(x)` or `pdf(x)` when weighted evaluations are needed.
+
+For SciPy-style multivariate normal objects, a lightweight adapter maps low-discrepancy points through `norm.ppf` and then injects covariance using a Cholesky factorization. This provides a stable, vectorized path for correlated Gaussian sampling.
+
+##### 2. Custom univariate mode with sanity checks
+
+Custom univariate objects are accepted when they follow a SciPy-like interface (`ppf` and optionally `pdf`/`logpdf`). At initialization, lightweight checks help identify common modeling mistakes:
+
+1. Non-monotone behavior in `ppf`.
+2. Non-finite values in `ppf`, `pdf`, or `logpdf`.
+3. Negative density values.
+4. Evidence that density behavior is inconsistent with normalization.
+
+These checks intentionally produce warnings rather than hard failures, preserving flexibility for advanced users while reducing silent numerical errors.
+
+#### Code Changes at a Glance
+
+The user-facing interface remains familiar.
+
+##### Existing independent workflow (unchanged)
+
+```python
+import scipy.stats as stats
+from qmcpy.discrete_distribution import DigitalNetB2
+from qmcpy.true_measure import SciPyWrapper
+
+tm = SciPyWrapper(
+ sampler=DigitalNetB2(2, seed=7),
+ scipy_distribs=[stats.norm(), stats.gamma(a=3.0)],
+)
+x = tm(4096)
+```
+
+##### New dependent-joint workflow
+
+```python
+import scipy.stats as stats
+from qmcpy.discrete_distribution import DigitalNetB2
+from qmcpy.true_measure import SciPyWrapper
+
+mvn = stats.multivariate_normal(
+ mean=[0.0, 0.0],
+ cov=[[1.0, 0.8], [0.8, 1.0]],
+)
+tm_joint = SciPyWrapper(DigitalNetB2(2, seed=7), scipy_distribs=mvn)
+x_joint = tm_joint(4096)
+```
+
+##### New custom-univariate workflow
+
+```python
+from qmcpy.discrete_distribution import DigitalNetB2
+from qmcpy.true_measure.triangular import TriangularDistribution
+from qmcpy.true_measure import SciPyWrapper
+
+tri = TriangularDistribution(c=0.3, loc=-1.0, scale=2.0)
+tm_custom = SciPyWrapper(DigitalNetB2(1, seed=11), scipy_distribs=tri)
+x_custom = tm_custom(4096)
+```
+
+#### Experimental Results
+
+The following six examples illustrate the new functionality in practice.
+
+##### Example 1: Independent vs dependent Gaussian behavior
+
+Independent normal marginals are compared against a dependent bivariate normal with target correlation `rho = 0.8`.
+
+* Figure 1: Independent vs dependent normals.
+
+
+
+The dependent construction reproduces both target correlation and mixed-moment behavior, while the independent baseline remains near zero correlation.
+
+##### Example 2: Zero-inflated exponential-uniform joint model
+
+A joint construction with an atom at `X=0` is used to model zero inflation with branch-dependent behavior in `Y`.
+
+* Figure 2: Zero-inflated exponential plus uniform joint model.
+
+
+
+The empirical mass at zero and the observed joint geometry are consistent with the model specification.
+
+##### Example 3: Acceptance-rejection dependence with MC vs QMC
+
+A triangular target region is sampled via acceptance-rejection using IID MC and digital-net QMC proposals.
+
+* Figure 3: Acceptance-rejection target with MC and QMC proposals.
+
+
+
+Both methods sample the intended target, while QMC provides visibly more even coverage and competitive moment estimates.
+
+##### Example 4: Custom triangular marginal
+
+A user-defined triangular distribution is passed directly to `SciPyWrapper`.
+
+* Figure 4: Custom triangular marginal with analytical overlay.
+
+
+
+Sample range and empirical behavior align closely with analytical expectations.
+
+##### Example 5: Warning path for malformed custom distributions
+
+An intentionally malformed custom distribution triggers diagnostic warnings, demonstrating that invalid definitions can be detected early rather than silently propagated.
+
+##### Example 6: Dependent multivariate Student t
+
+A dependent multivariate Student t model is used to validate heavy-tail behavior beyond the Gaussian setting.
+
+* Figure 5: Dependent multivariate Student t behavior.
+
+
+
+Empirical dependence and covariance trends are consistent with the target model and expected heavy-tailed structure.
+
+#### Validation
+
+The behavior above is supported by dedicated unit checks, including:
+
+1. Correlation and mixed moments in joint Gaussian sampling.
+2. Zero-mass accuracy in zero-inflated sampling.
+3. Range and mean behavior for custom triangular marginals.
+4. Correlation preservation in dependent Student t sampling.
+
+#### Practical Impact
+
+This extension keeps the familiar QMCPy workflow while enabling substantially richer probabilistic modeling. Users can now encode dependence, prototype domain-specific marginals, and retain diagnostic guardrails inside the same interface that previously handled only independent SciPy marginals.
+
+#### Conclusion
+
+The updated `SciPyWrapper` preserves simplicity for existing users and adds the expressive power needed for realistic modern simulation tasks. By supporting dependent joints, custom marginals, and warning-based validation, it strengthens both modeling flexibility and numerical reliability within QMCPy's core workflow.
+
+#### References
+
+1. Choi, S.-C., Hickernell, F., McCourt, M., Sorokin, A. QMCPy: A Quasi-Monte Carlo Python Library. https://qmcsoftware.github.io/QMCSoftware/
+2. SciPy Statistics Documentation. https://docs.scipy.org/doc/scipy/reference/stats.html
+3. QMCSoftware repository. https://github.com/QMCSoftware/QMCSoftware
diff --git a/blog/speeding-up-qmcpy-with-distributable-c-code/index.qmd b/blog/speeding-up-qmcpy-with-distributable-c-code/index.qmd
new file mode 100644
index 0000000..8c23393
--- /dev/null
+++ b/blog/speeding-up-qmcpy-with-distributable-c-code/index.qmd
@@ -0,0 +1,124 @@
+---
+title: "Speeding up QMCPy with Distributable C Code"
+author: "Aleksei Sorokin and Jagadeeswaran Rathinavel"
+date: 2021-03-21
+date-format: "MMMM D, YYYY"
+description: "Using distributable C extensions through ctypes and packaging metadata to accelerate QMCPy."
+categories:
+ - "QMCPy"
+ - "Performance"
+ - "C Extensions"
+---
+
+This post explains how QMCPy uses distributable C extensions through `ctypes` and packaging metadata to speed up low-discrepancy generators.
+
+Many Python packages rely on underlying C or C++ code to speed up their numerical methods. For example, [NumPy](https://numpy.org/) calls C and C++ extensions in order to speed up matrix manipulation algorithms.
+
+Real Python's article [*Python Bindings: Calling C or C++ From Python*](https://realpython.com/python-bindings-overview/#python-bindings-overview) discusses a few reasons why you may want to utilize C or C++ extensions within your Python package. Perhaps you already have a stable library in C or C++ that you want to call from Python. Our approach in this blog will allow the existing extension to be called from Python with only minor code modifications. You may also be interested in speeding up your Python code by moving it to a compiled language that can optimize subroutines. For example, in QMCPy we have found that many low discrepancy sequence generators are significantly faster when implemented in C.
+
+So why not implement everything in C or C++? In our experience, Python delivers a convenience, readability, and community engagement that allow for rapid development, testing, and distribution to a large audience of active users.
+
+While the benefits of moving certain modules to C/C++ have been well documented, implementing these extensions to play nicely with your existing Python codebase can often be quite tricky. Moreover, writing extensions for platform-independent distribution with PyPI, so someone can `pip install yourPackage`, can be even more challenging. In this blog post we share how we developed the QMCPy package [1] to be platform agnostic while utilizing C extensions.
+
+## C
+
+When you first explore writing C/C++ extensions you will likely come across [Python's recommended solution](https://docs.python.org/3/extending/extending.html). This approach requires a good bit of boilerplate code and prohibits plug-and-play of an existing C/C++ library. QMCPy's approach uses the [`ctypes`](https://docs.python.org/3/library/ctypes.html) library to call a C function with a few lines of Python defining the arguments and return values of the compiled function.
+
+Let us now turn to an example from our QMCPy library. Based on Art Owen's work in [2], we wrote the below Halton generator in C. Note that the implementation does not contain all the boilerplate code of a native Python solution, but instead may be used as a standalone C file.
+
+```c
+#include "MRG63k3a.h"
+static int primes[1000] = {2, 3, 5, 7, ...};
+EXPORT void halton_owen(int n, int d, int n0, int d0,
+ int randomize, double *ans, long long seed){
+ seed_MRG63k3a(seed); ...}
+```
+
+A few important notes about the above code are the use of `#include`, `EXPORT`, and `long long`. Depending on the compiler, such as `gcc` or Windows `cl.exe`, `EXPORT` allows us to expose a function, in the above case `halton_owen`, so that the Python code can invoke it. When you `EXPORT` a function it makes the C code available to `ctypes`. The Halton generator utilizes the MRG63k3a random number generator [3], which is stored in a separate file. We can call this function by creating a `.h` file that defines the external function we wish to call. In `MRG63k3a.h` we define the `seed_MRG63k3a` method which is then included and used in the above Halton generator.
+
+When your Python package is installed, the C compiler that builds the extensions is platform-specific. We found that the `gcc` compiler uses 8 bytes to store a `long` while Windows `cl.exe` uses only 4. As a workaround, we suggest using the `long long` datatype, which is 8 bytes for both `gcc` and `cl.exe`. A nice way to verify you are using `gcc` and debug these cross-language problems is to intentionally trigger compiler errors.
+
+With these three C files, the Halton generator, MRG63k3a, and MRG63k3a's header, we are ready to call our function from Python.
+
+## Python Code
+
+First, we will use `ctypes` to define our function from Python. `ctypes` requires that we define the arguments and return values of our Halton function in order for it to be treated like a native Python function. Below is an example of how to set up and call our Halton function in C using Python.
+
+```python
+import ctypes
+from ctypes import CDLL, RTLD_GLOBAL
+from os.path import dirname, abspath
+from glob import glob
+from numpy import *
+
+# load the library
+c_lib = CDLL(
+ glob(dirname(abspath(__file__))+'/c_lib*')[0],
+ mode = RTLD_GLOBAL)
+
+# define the function arguments
+halton_cf = c_lib.halton_owen
+halton_cf.argtypes = [
+ ctypes.c_int, # n
+ ctypes.c_int, # d
+ ctypes.c_int, # n0
+ ctypes.c_int, # d0
+ ctypes.c_int, # randomize
+ ctypeslib.ndpointer( # result array
+ ctypes.c_double,
+ flags = 'C_CONTIGUOUS'),
+ ctypes.c_long] # seed
+
+# define the return value
+halton_cf.restype = None
+
+# example call to the function
+# create an empty array to fill with Halton points
+x = zeros((5,3), dtype=double)
+# fill the array with 5, 3-dimensional Halton points
+halton_cf(5, 3, 0, 0, True, x, 17)
+```
+
+The second piece of Python code you will need is a `setup.py`. The `setup.py` file defines the C extensions of your package and helps prepare your package for distribution on PyPI. While it is possible to compile and call your extension function without a `setup.py` file, we found this method to be the easiest and most straightforward for package distribution.
+
+Below is a snippet from our `setup.py` file that defines the extensions, packages, and other metadata. Note that we use the [`setuptools`](https://setuptools.readthedocs.io/en/latest/) package to easily define our distribution properties, although [`distutils`](https://docs.python.org/3/library/distutils.html) may also be used.
+
+```python
+import setuptools
+from setuptools import Extension
+
+# define package API
+packages = [
+ 'qmcpy',
+ 'qmcpy.discrete_distribution',
+ 'qmcpy.discrete_distribution.halton']
+
+setuptools.setup(
+ name="qmcpy",
+ packages=packages,
+ include_package_data=True,
+ ext_modules=[
+ Extension(
+ name='qmcpy.discrete_distribution.c_lib.c_lib',
+ sources=[
+ 'qmcpy/discrete_distribution/c_lib/halton_owen.c',
+ 'qmcpy/discrete_distribution/c_lib/MRG63k3a.c'],
+ )],)
+```
+
+When distributing your package on PyPI, you may come across an error regarding missing `.h` files, e.g., the `MRG63k3a.h` file mentioned earlier. For this file to be included in your distribution, you need to include a `MANIFEST.in` file that defines all the non-Python and non-C code to be included. That way the `.h` files and other files will be included in your package distribution. We provide a sample from our `MANIFEST.in` below.
+
+```text
+include qmcpy/discrete_distribution/c_lib/*.h
+include qmcpy/discrete_distribution/sobol/generating_matrices/*.npy
+include qmcpy/discrete_distribution/lattice/generating_vectors/*.npy
+```
+
+## References
+
+1. Choi, S.-C. T., Hickernell, F., McCourt, M., & Sorokin, A. QMCPy: A
+quasi-Monte Carlo Python Library. [https://qmcsoftware.github.io/QMCSoftware/](https://qmcsoftware.github.io/QMCSoftware/). 2020.
+2. Owen, A. B. A randomized Halton algorithm in R. 2017.
+ [arXiv:1706.02808 [stat.CO]](https://arxiv.org/abs/1706.02808).
+3. L'Ecuyer, P. Good parameters and implementations for combined
+multiple recursive random number generators. *Operations Research*, 47, 159-164. [https://pubsonline.informs.org/doi/abs/10.1287/opre.47.1.159](https://pubsonline.informs.org/doi/abs/10.1287/opre.47.1.159). 1999.
diff --git a/blog/visualizing-the-generated-samples-helps/figures/halton-uniform.png b/blog/visualizing-the-generated-samples-helps/figures/halton-uniform.png
new file mode 100644
index 0000000..fd3cfcf
Binary files /dev/null and b/blog/visualizing-the-generated-samples-helps/figures/halton-uniform.png differ
diff --git a/blog/visualizing-the-generated-samples-helps/figures/halton.png b/blog/visualizing-the-generated-samples-helps/figures/halton.png
new file mode 100644
index 0000000..ebc87c5
Binary files /dev/null and b/blog/visualizing-the-generated-samples-helps/figures/halton.png differ
diff --git a/blog/visualizing-the-generated-samples-helps/figures/iid-gaussian.png b/blog/visualizing-the-generated-samples-helps/figures/iid-gaussian.png
new file mode 100644
index 0000000..f74f47d
Binary files /dev/null and b/blog/visualizing-the-generated-samples-helps/figures/iid-gaussian.png differ
diff --git a/blog/visualizing-the-generated-samples-helps/figures/iid.png b/blog/visualizing-the-generated-samples-helps/figures/iid.png
new file mode 100644
index 0000000..64a3379
Binary files /dev/null and b/blog/visualizing-the-generated-samples-helps/figures/iid.png differ
diff --git a/blog/visualizing-the-generated-samples-helps/index.qmd b/blog/visualizing-the-generated-samples-helps/index.qmd
new file mode 100644
index 0000000..ef00069
--- /dev/null
+++ b/blog/visualizing-the-generated-samples-helps/index.qmd
@@ -0,0 +1,102 @@
+---
+title: "Visualizing the Generated Samples Helps"
+author: "Aadit Jain"
+date: 2026-02-25
+date-format: "MMMM D, YYYY"
+description: "Using QMCPy's plot_proj function to visualize discrete distributions and true measures."
+categories:
+ - "QMCPy"
+ - "Visualization"
+ - "Sampling"
+image: figures/halton-uniform.png
+---
+
+This post introduces QMCPy's `plot_proj` function for visualizing two-dimensional projections of discrete distributions and true measures.
+
+It is a universal truth that visuals better appeal to the human mind than a group of numbers listed out. With visuals, it is easier for us to discern patterns and identify any flaws in the logic behind our mathematical calculations and equations. To visualize the different Discrete Distribution and True Measure objects, the [QMCPy](https://qmcpy.org/) Plot Projections function has been developed. This blog presents the different applications of this function.
+
+## What Does the Plot Projections Function Do?
+
+A Discrete Distribution or True Measure object with $d \geq 2$
+dimensions has a maximum of $d\times(d -1)$ dimensional pairings. For
+example, $[2,3]$ and $[3,2]$ are considered separate pairings. The
+Plot Projections function, referred to as `plot_proj`, plots all or a
+subset of all the possible dimension pairings based on user arguments.
+This function either takes a Discrete Distribution or True Measure object
+at a time. It can also display extensibility by passing in a list of
+successively larger samples, for example $[2^6,2^7,2^8]$.
+
+To display extensibility, the `plot_proj` function uses [the default `prop_cycle`](https://matplotlib.org/stable/gallery/color/color_cycle_default.html), which is obtained from the `rc` parameters of Matplotlib. This default `prop_cycle` contains a list of colors through which the `plot_proj` function iterates over and displays extensibility. The list of colors is: blue, orange, green, red, purple, brown, pink, grey, yellow, cyan. The colors are stored in this order but in a hexadecimal format.
+
+The parameters and plot examples of this function can be seen in the [Plot Projections Notebook](https://github.com/QMCSoftware/QMCSoftware/blob/master/demos/plot_proj_function.ipynb).
+
+## Setting Up the QMCPy Environment Before Utilizing the Plot Projections Function
+
+```python
+import qmcpy as qp
+```
+
+## The Different Applications of the Plot Projection Function
+
+1. Here we show a two dimensional projection of an IID object:
+
+ ```python
+ d = 2
+ iid = qp.IIDStdUniform(d)
+ fig,ax = qp.plot_proj(iid, n = 2**7)
+ ```
+
+ **Figure 1: Uniform IID Object Projection**
+
+ 
+
+2. Here we show a two dimensional projection of a Gaussian object and how
+the axes returned by the `plot_proj` function can be manipulated by adding a horizontal and vertical line to denote the x and y axis respectively:
+
+ ```python
+ d = 2
+ iid = qp.IIDStdUniform(d)
+ iid_gaussian = qp.Gaussian(iid,mean =[2,4],covariance=[[9,4],[4,5]])
+ fig,ax = qp.plot_proj(iid_gaussian, n = 2**7)
+ ax[0,0].axvline(x=0,color= 'k',alpha=.25); #adding vertical line
+ ax[0,0].axhline(y=0,color= 'k',alpha=.25); #adding horizontal line
+ ```
+
+ **Figure 2: Gaussian IID Object Projection**
+
+ 
+
+3. Here we show certain specified dimensional projections, with
+dimensions 1 and 2 on the x axes and dimensions 3 and 4 on the y axes, of a Uniform object with successively increasing numbers of points. The initial points are in blue. The next additional points are in orange. The final additional points are in green:
+
+ ```python
+ d = 4
+ halton = qp.Halton(d)
+ halton_uniform = qp.Uniform(halton,lower_bound=[1,2,3,4],
+ upper_bound=[5,7,9,11])
+ fig, ax = qp.plot_proj(halton_uniform, n = [2**6, 2**7, 2**8],
+ d_horizontal = [1,2], d_vertical = [3,4])
+ ```
+
+ **Figure 3: Halton Object Projection**
+
+ 
+
+4. Here we show a four dimensional projection of a Halton object with
+successively increasing numbers of points. The initial points are in blue. The next additional points are in orange. The final additional points are in green:
+
+ ```python
+ d = 4
+ halton = qp.Halton(d)
+ fig,ax = qp.plot_proj(halton, n = [2**5, 2**6, 2**7],
+ d_horizontal = range(d), d_vertical = range(d),
+ math_ind = False, marker_size = 15)
+ ```
+
+ **Figure 4: Halton Object Projection (More Points)**
+
+ 
+
+## How This Function Benefits Us
+
+In addition to making it easy to see the difference between different Discrete Distribution and True Measure objects, this function consists of many features that makes it user-friendly and help generate a strong and precise visualization of the different Discrete Distribution and True Measure objects. For instance, the extensibility feature enables us to see how the space of the plot fills up, the `marker_size` parameter helps make the samples/points bigger when plotting high-dimensional projections, and the `math_ind` parameter allows the user to either input mathematical or Python dimensions for the sampler based on one's preference. This function could also be developed in the future to support other distributions such as Brownian Motion.
diff --git a/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/accumulate_data_blog_uml.png b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/accumulate_data_blog_uml.png
new file mode 100644
index 0000000..48beb4b
Binary files /dev/null and b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/accumulate_data_blog_uml.png differ
diff --git a/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/discrete_distribution_blog_uml.png b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/discrete_distribution_blog_uml.png
new file mode 100644
index 0000000..230208a
Binary files /dev/null and b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/discrete_distribution_blog_uml.png differ
diff --git a/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/integrand_blog_uml.png b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/integrand_blog_uml.png
new file mode 100644
index 0000000..d824ee2
Binary files /dev/null and b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/integrand_blog_uml.png differ
diff --git a/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/qmcpy_blog_uml3.png b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/qmcpy_blog_uml3.png
new file mode 100644
index 0000000..7e6e3f8
Binary files /dev/null and b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/qmcpy_blog_uml3.png differ
diff --git a/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/qmcpy_uml1-1.png b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/qmcpy_uml1-1.png
new file mode 100644
index 0000000..521c8a9
Binary files /dev/null and b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/qmcpy_uml1-1.png differ
diff --git a/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/qmcpy_uml2.png b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/qmcpy_uml2.png
new file mode 100644
index 0000000..797be69
Binary files /dev/null and b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/qmcpy_uml2.png differ
diff --git a/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/stopping_criterion_blog_uml.png b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/stopping_criterion_blog_uml.png
new file mode 100644
index 0000000..55e3e21
Binary files /dev/null and b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/stopping_criterion_blog_uml.png differ
diff --git a/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/true_measure_uml.png b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/true_measure_uml.png
new file mode 100644
index 0000000..8cd280b
Binary files /dev/null and b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/true_measure_uml.png differ
diff --git a/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/util_err.png b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/util_err.png
new file mode 100644
index 0000000..fa39c55
Binary files /dev/null and b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/util_err.png differ
diff --git a/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/util_warn.png b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/util_warn.png
new file mode 100644
index 0000000..f35c43f
Binary files /dev/null and b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/figures/util_warn.png differ
diff --git a/blog/visualizing-the-internals-of-object-classes-in-qmcpy/index.qmd b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/index.qmd
new file mode 100644
index 0000000..5328933
--- /dev/null
+++ b/blog/visualizing-the-internals-of-object-classes-in-qmcpy/index.qmd
@@ -0,0 +1,118 @@
+---
+title: "Visualizing the Internals of Object Classes in QMCPy"
+author: "Sou-Cheng Choi and Aleksei Sorokin"
+date: 2021-02-25
+date-format: "MMMM D, YYYY"
+description: "UML diagrams explaining QMCPy's object-oriented architecture and core class relationships."
+categories:
+ - "QMCPy"
+ - "Software Architecture"
+ - "UML"
+image: figures/qmcpy_uml1-1.png
+---
+
+This post uses UML diagrams to explain QMCPy's object-oriented architecture and relationships among its core classes.
+
+As a software library grows, so does its complexity. This comment certainly applies to QMCPy [1], our Python library for high-dimensional numerical integration. [UML (Unified Modelling Language) diagrams](https://en.wikipedia.org/wiki/Unified_Modeling_Language) are a helpful tool for visualizing QMCPy's intricate object-oriented framework. These network diagrams display an object's methods, attributes, dependencies, and inheritance relationships. We have used the Python tool [`pyreverse`](https://pypi.org/project/pyreverse/) to automatically generate such UML diagrams, which we have included in the [QMCPy documentation](https://qmcsoftware.github.io/QMCSoftware/). For a comprehensive introduction to various aspects of UML and its latest version 2.5, readers may refer to, for example, [2].
+
+## Overview of QMCPy Classes
+
+First, we overview the relationships between the five main abstract classes in [QMCPy version 1.0](https://qmcpy.org/2021/02/12/qmcpy-version-1-0/):
+
+- `Integrand`,
+- `TrueMeasure`,
+- `DiscreteDistribution`,
+- `StoppingCriterion`, and
+- `AccumulateData`.
+
+For clearer illustration and better readability, we may not include all subclasses implemented in QMCPy in the subsequent diagrams. Interested readers are referred again to the [QMCPy documentation](https://qmcsoftware.github.io/QMCSoftware/).
+
+In a UML class diagram, each class is contained in a rectangular box. A child class has an edge with a triangular arrow head that points to its parent. If a class **C** internally uses an object of another class **D**, the edge would have a solid black diamond head from class **D** pointing to **C**. A green label of an edge recaps the name of a field in the class being pointed to, realized by the class from which the edge stems.
+
+The first UML class diagram shows the abstract `Integrand` class and its implementations (children). Notice that many integrands used to price financial options, specifically `AsianOption`, `EuropeanOption`, and `MLCallOptions`, utilize `BrownianMotion`, a child of the `Gaussian` class and grandchild of the abstract `TrueMeasure` class. In particular, `AsianOption` and `EuropeanOption` contain a field called `true_measure`, which is highlighted in green in the UML class diagram below and implemented as a `BrownianMotion` object.
+
+
+
+ Overview UML class diagram for the abstract Integrand class and selected implementations.
+
+
+The second UML class diagram is `DiscreteDistribution` and its subclasses.
+
+
+
+ Overview UML class diagram for DiscreteDistribution and its subclasses.
+
+
+The last high-level diagram relates `StoppingCriterion` and `AccumulateData`. In particular, every `StoppingCriterion` implementation uses an `AccumulateData` implementation for storing the parameters that were used and set during the numerical approximation algorithm. The following class diagram includes only Quasi-Monte Carlo stopping criteria, but QMCPy actually also contains a number of standard (IID) Monte Carlo stopping criteria as well.
+
+
+
+ Overview UML class diagram relating StoppingCriterion and AccumulateData.
+
+
+## More Details of QMCPy Classes
+
+In the remainder of this blog, we will present in greater detail the internal members of each main class. Each class is listed at the top of a rectangular box with its public fields and methods in the middle and bottom sections of the box, respectively. A child class inherits the methods of its parent class. However, a child class may override the parent's handed down method. In this case, the child class method is listed again at the bottom of its UML box.
+
+The `Integrand` class has three main fields and methods. Each of its
+five subclasses has its own specific implementation of the integrand,
+$g$, but shares the same method, $f$, that returns the weighted average
+of the transformed function at nodes from the `DiscreteDistribution`.
+The specific weights, transformation, and sampling mechanism are
+determined by the associated realizations of the `TrueMeasure` and
+`DiscreteDistribution` classes. For instance, the `Keister` integrand
+uses the `Gaussian` true measure and may be paired with any
+`DiscreteDistribution` instance.
+
+
+
+ Detailed UML class diagram for Integrand.
+
+
+QMCPy has implemented five children classes for `TrueMeasure`. A child class of `TrueMeasure` has an attribute called `discrete_distrib` that is a `DiscreteDistribution` instance. This enables the main `gen_samples` method to select and transform points accordingly.
+
+
+
+ Detailed UML class diagram for TrueMeasure.
+
+
+`DiscreteDistribution` in QMCPy plays a central role in (Q)MC algorithms, which are iterative in nature. In every iteration, a concrete subclass of `DiscreteDistribution` decides the coordinates of sampling points for integrand evaluations, which are then aggregated into an average value that serves as an estimate of a given integral problem. We refer readers to an earlier blog for a succinct presentation of [low discrepancy sampling points](../what-makes-a-sequence-low-discrepancy/index.qmd) used in QMC algorithms versus IID sampling schemes in more traditional MC methods.
+
+
+
+ Detailed UML class diagram for DiscreteDistribution.
+
+
+QMCPy's abstract `StoppingCriterion` class currently has the largest number of instances. Each concrete implementation has two main abstract methods, `set_tolerance` and `integrate`. The method `set_tolerance` allows users to reset absolute and/or relative tolerances used in the `integrate` method. Calling `integrate` will construct an `AccumulateData` object to generate and house data such as sampling indices, function evaluations, and expectation approximations.
+
+
+
+ Detailed UML class diagram for StoppingCriterion.
+
+
+As mentioned, an `AccumulateData` subclass collects data throughout the numerical integration computation. The method `update_data` collects statistical estimates such as the sample mean, sample variance, and approximate time per sample. An `AccumulateData` instance can often be used by multiple `StoppingCriterion`. For example, `LDTransformData` is used by both the `CubQMCLatticeG` and `CubQMCSobolG` stopping criteria. Since an `AccumulateData` object knows about the four other components in the QMC problem, printing the data object displays a nice summary of all relevant fields and parameters.
+
+
+
+ Detailed UML class diagram for AccumulateData.
+
+
+Lastly, QMCPy has extended Python's `Warning` and `Exception` classes to provide developers specific types of errors and warnings.
+
+
+
+ UML class diagram for QMCPy-specific exceptions.
+
+
+
+
+ UML class diagram for QMCPy-specific warnings.
+
+
+We hope these UML diagrams help both researchers and developers better understand the QMCPy architecture. The diagrams throughout this blog are reproducible using the `pyreverse` package whose `S` option will reveal extra details including private fields and methods.
+
+## References
+
+1. Choi, S.-C. T., Hickernell, F., McCourt, M., & Sorokin, A. QMCPy: A
+quasi-Monte Carlo Python Library. [https://qmcsoftware.github.io/QMCSoftware/](https://qmcsoftware.github.io/QMCSoftware/). 2020.
+2. Unhelkar, B. *Software Engineering with UML*. CRC Press, 2017.
diff --git a/blog/what-makes-a-sequence-low-discrepancy/figures/iid_uniform_pts.png b/blog/what-makes-a-sequence-low-discrepancy/figures/iid_uniform_pts.png
new file mode 100644
index 0000000..c82f010
Binary files /dev/null and b/blog/what-makes-a-sequence-low-discrepancy/figures/iid_uniform_pts.png differ
diff --git a/blog/what-makes-a-sequence-low-discrepancy/figures/lattice_pts.png b/blog/what-makes-a-sequence-low-discrepancy/figures/lattice_pts.png
new file mode 100644
index 0000000..0bb525e
Binary files /dev/null and b/blog/what-makes-a-sequence-low-discrepancy/figures/lattice_pts.png differ
diff --git a/blog/what-makes-a-sequence-low-discrepancy/index.qmd b/blog/what-makes-a-sequence-low-discrepancy/index.qmd
new file mode 100644
index 0000000..1d01683
--- /dev/null
+++ b/blog/what-makes-a-sequence-low-discrepancy/index.qmd
@@ -0,0 +1,219 @@
+---
+title: "What Makes a Sequence Low Discrepancy?"
+author: "Fred J. Hickernell"
+date: 2020-07-08
+date-format: "MMMM D, YYYY"
+description: "Discrepancy as a measure of uniformity and why low-discrepancy sequences improve QMC integration."
+categories:
+ - "Quasi-Monte Carlo"
+ - "Low-Discrepancy Sequences"
+ - "Discrepancy"
+image: figures/lattice_pts.png
+---
+
+This post introduces discrepancy as a way to measure uniformity and explains why low-discrepancy sequences improve QMC integration.
+
+The first blog post, [Why Add Q to MC?](../why-add-q-to-mc/index.qmd), introduced the concept of evenly spread points, which are commonly referred to as *low discrepancy* (LD) points. This is in contrast to independent and identically distributed (IID) points.
+
+Consider two sequences,
+
+$$
+\boldsymbol{T}_1, \boldsymbol{T}_2, \ldots
+\overset{\text{IID}}{\sim} \mathcal{U}[0,1]^d,
+$$
+
+
+
+ Figure 1: 64 IID standard uniform points in 2 dimensions.
+
+
+and
+
+$$
+\boldsymbol{X}_1, \boldsymbol{X}_2, \ldots
+\overset{\text{LD}}{\sim} \mathcal{U}[0,1]^d.
+$$
+
+
+
+ Figure 2: 64 shifted lattice points in 2 dimensions.
+
+
+Both sequences are expected to look like points spread uniformly over
+the unit cube, $[0,1]^d$. The first sequence must be random, or as
+random looking as our deterministic random number generators can make
+it. Since the points are independent, the location of any
+$\boldsymbol{T}_i$ has no bearing on the location of any other
+$\boldsymbol{T}_j$. Removing a point at random does not affect the IID
+property.
+
+The second sequence may be random or deterministic. Let
+$F_{\{\boldsymbol{X}_1, \ldots, \boldsymbol{X}_n\}}$ denote the
+empirical distribution function of the first $n$ points of this
+sequence, i.e., the probability distribution that assigns a probability
+of $1/n$ to each location $\boldsymbol{X}_i$. For
+$\boldsymbol{X}_1, \boldsymbol{X}_2, \ldots$ to be LD,
+$F_{\{\boldsymbol{X}_1, \ldots, \boldsymbol{X}_n\}}$ should be close to
+the uniform probability distribution,
+$F_{\text{unif}}: \boldsymbol{x} \mapsto x_1 \cdots x_d$.
+
+"Close" implies that we can measure how far apart two distributions are.
+We call this measure the *discrepancy*. Just like beauty is in the eye of
+the beholder, there are different measurements of discrepancy. They tend
+to take the form of the distance between the empirical distribution of
+the point set, $F_{\{\boldsymbol{X}_1, \ldots, \boldsymbol{X}_n\}}$, and
+the target measure, $F_{\text{unif}}$. An example is the star
+discrepancy [1, (3.16)]:
+
+$$
+\operatorname{disc}(\{\boldsymbol{X}_1, \ldots, \boldsymbol{X}_n\})
+:=
+\sup_{\boldsymbol{x} \in [0,1]^d}
+\left|
+F_{\text{unif}}(\boldsymbol{x})
+- F_{\{\boldsymbol{X}_1, \ldots, \boldsymbol{X}_n\}}(\boldsymbol{x})
+\right|.
+$$
+
+This quantity is known in the statistics literature as a
+Kolmogorov-Smirnov goodness-of-fit statistic. This discrepancy is the
+maximum absolute difference between the volume of the box
+$[\boldsymbol{0}, \boldsymbol{x}]^d$ and the proportion of the points
+$\{\boldsymbol{X}_1, \ldots, \boldsymbol{X}_n\}$ that lie in that box.
+Ideally, these should be the same, but practically they will be at least
+a bit different.
+
+The computational cost of evaluating the star discrepancy can be rather
+large, typically at least $\mathcal{O}(n^d)$ operations. A family of
+computationally cheaper discrepancies is defined in terms of a kernel,
+$K: [0,1]^d \times [0,1]^d \to \mathbb{R}$, which satisfies two crucial
+properties:
+
+$$
+\begin{aligned}
+\text{Symmetry:} \quad&
+K(\boldsymbol{t}, \boldsymbol{x})
+= K(\boldsymbol{x}, \boldsymbol{t})
+\qquad \forall \boldsymbol{t}, \boldsymbol{x} \in [0,1]^d, \\
+\text{Positive definiteness:} \quad&
+\boldsymbol{c}^{\mathsf{T}} \mathsf{K} \boldsymbol{c} > 0,
+\text{ where }
+\mathsf{K}
+= \bigl(K(\boldsymbol{x}_i, \boldsymbol{x}_j)\bigr)_{i,j=1}^n, \\
+& \qquad \qquad
+\forall \boldsymbol{c} \ne \boldsymbol{0},
+\text{ distinct } \boldsymbol{x}_1, \boldsymbol{x}_2, \ldots
+\in [0,1]^d.
+\end{aligned}
+$$
+
+For such a kernel, we may define a discrepancy as
+
+$$
+\begin{aligned}
+&\operatorname{disc}(\{\boldsymbol{X}_1, \ldots, \boldsymbol{X}_n\}) \\
+&\quad :=
+\int_{[0,1]^d \times [0,1]^d}
+K(\boldsymbol{t}, \boldsymbol{x})
+\, \mathrm{d}
+\left(F_{\text{unif}}
+- F_{\{\boldsymbol{X}_1, \ldots, \boldsymbol{X}_n\}}\right)
+(\boldsymbol{t})
+\, \mathrm{d}
+\left(F_{\text{unif}}
+- F_{\{\boldsymbol{X}_1, \ldots, \boldsymbol{X}_n\}}\right)
+(\boldsymbol{x}) \\
+&\quad =
+\int_{[0,1]^d \times [0,1]^d}
+K(\boldsymbol{t}, \boldsymbol{x})
+\, \mathrm{d}\boldsymbol{t}\,\mathrm{d}\boldsymbol{x}
+- \frac{2}{n}
+\sum_{i=1}^n
+\int_{[0,1]^d}
+K(\boldsymbol{x}_i, \boldsymbol{x})
+\, \mathrm{d}\boldsymbol{x} \\
+&\qquad
++ \frac{1}{n^2}
+\sum_{i,j=1}^n K(\boldsymbol{x}_i, \boldsymbol{x}_j).
+\end{aligned}
+$$
+
+For example, the centered $L^2$-discrepancy [2] is defined in terms of
+the kernel
+
+$$
+K(\boldsymbol{t}, \boldsymbol{x})
+=
+\prod_{k=1}^d
+\left[
+1
++ \frac{1}{2}|t_k - 1/2|
++ \frac{1}{2}|x_k - 1/2|
+- \frac{1}{2}|t_k - x_k|
+\right].
+$$
+
+After straightforward calculations it becomes
+
+$$
+\begin{aligned}
+&\operatorname{disc}(\{\boldsymbol{X}_1, \ldots, \boldsymbol{X}_n\})
+=
+\left(\frac{13}{12}\right)^d
+- \frac{2}{n}
+\sum_{i=1}^n
+\prod_{k=1}^d
+\left(
+1
++ \frac{1}{2}|x_{ik} - 1/2|
+- \frac{1}{2}|x_{ik} - 1/2|^2
+\right) \\
+&\qquad
++ \frac{1}{n^2}
+\sum_{i,j=1}^n
+\prod_{k=1}^d
+\left[
+1
++ \frac{1}{2}|x_{ik} - 1/2|
++ \frac{1}{2}|x_{jk} - 1/2|
+- \frac{1}{2}|x_{ik} - x_{jk}|
+\right].
+\end{aligned}
+$$
+
+This discrepancy only requires $\mathcal{O}(dn^2)$ operations to
+evaluate.
+
+LD sequences have discrepancies of $\mathcal{O}(n^{-1+\epsilon})$ for
+the discrepancies illustrated above. IID sequences have root mean square
+discrepancies of $\mathcal{O}(n^{-1/2})$. This difference in asymptotic
+order can translate into orders of magnitude improvements in the
+accuracy of numerical solutions.
+
+For problems where $d$ is large, the discrepancies defined above do not
+decay so quickly. However, if these discrepancy definitions are modified
+to include *coordinate weights* [1, Section 4], then they retain their
+$\mathcal{O}(n^{-1+\epsilon})$ decay. Coordinate weights express the
+assumption that certain coordinates contribute more to the variation of
+the function than others.
+
+Demonstrating that a particular sequence is LD can be done by brute
+force computation, which requires in general $\mathcal{O}(dn^2)$
+operations. For certain sequences matched with certain discrepancy
+definitions, this can be reduced to $\mathcal{O}(dn)$ operations. If $n$
+is small enough, the search for an LD set can be performed using global
+optimization algorithms [3]. Number theoretic arguments are used to
+construct certain popular LD sequences [4, 5].
+
+## References
+
+1. Dick, J., Kuo, F., & Sloan, I. H. High dimensional integration: The
+ quasi-Monte Carlo way. *Acta Numerica*, 22, 133-288 (2013).
+2. Hickernell, F. J. A generalized discrepancy and quadrature error
+ bound. *Mathematics of Computation*, 67, 299-322 (1998).
+3. Winker, P., & Fang, K. T. Application of threshold accepting to the
+evaluation of the discrepancy of a set of points. *SIAM Journal on Numerical Analysis*, 34, 2028-2042 (1997).
+4. Dick, J., & Pillichshammer, F. *Digital Nets and Sequences:
+Discrepancy Theory and Quasi-Monte Carlo Integration*. Cambridge University Press, Cambridge (2010).
+5. Niederreiter, H. *Random Number Generation and Quasi-Monte Carlo
+ Methods*. SIAM, Philadelphia (1992).
diff --git a/styles/custom.scss b/styles/custom.scss
index 12830e3..b9f1eb1 100644
--- a/styles/custom.scss
+++ b/styles/custom.scss
@@ -31,9 +31,33 @@ $border-radius: 0.4rem;
}
main.content {
+ min-width: 0;
+ overflow-x: clip;
font-size: 1.375rem;
}
+main.content section,
+main.content div.sourceCode,
+main.content pre.sourceCode {
+ min-width: 0;
+ max-width: 100%;
+}
+
+main.content div.sourceCode,
+main.content pre.sourceCode {
+ overflow-x: auto;
+}
+
+main.content figure img {
+ max-width: 100%;
+ height: auto;
+}
+
+main.content figure > div[style*="display: grid"] > img {
+ width: 100%;
+ min-width: 0;
+}
+
.hero {
margin: 0 0 2.5rem;
padding: clamp(2rem, 6vw, 4.5rem);
+ 
+ 
+ 
+ 
+ 
+
+
+ 
+ 
+ 
+
+
+
+
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+