VESTA: Visual Exploration with Statistical Tool Agents

# Vesta: Visual Exploration with Statistical Tool Agents
Source: [https://arxiv.org/html/2606.00384](https://arxiv.org/html/2606.00384)
William Rudman1, Abhishek Divekar1, Kanishk Jain1, Sebastian Joseph1, Stella S\. R\. Offner1, Matthew Lease1, Kyle Mahowald1, Greg Durrett2, Junyi Jessy Li1 1The University of Texas at Austin2New York University Correspondence:william\.rudman@utexas\.edu\|[https://github\.com/wrudman/VESTA](https://github.com/wrudman/VESTA)

###### Abstract

Fitting quantitative models to data is a central step in scientific workflows, yet it remains one of the least automated\. Recent agent\-based systems leverage language and vision\-language models \(VLMs\) to iteratively propose and refine statistical models, but these systems struggle on more challenging modeling tasks\. To address these limitations, we introduceVesta\(VisualExploration withStatisticalToolAgents\), a framework that equips VLMs with a dynamically growing exploration toolkit to guide model refinement through data transformations, hypothesis\-driven visualizations, and robust statistical tests\. Unlike prior systems that rely on iterative critique alone,Vestaactively explores data before and during refinement by selecting or creating diagnostic tools, which accumulate in the model’s context and can be reused later\. We evaluateVestaagainst established baselines in three toolkit configurations: no tools, static expert\-written tools, and dynamic model\-written tools\. To support this evaluation, we introduceDawn\(Dataset forAutomatedWorkflows andNumerical Modeling\), a benchmark targeting distribution fitting and time series modeling with varying difficulty tiers, and culminating in real\-world astronomy tasks including modeling initial mass functions and gravitational\-wave chirp signals\. We find thatVesta’s dynamic tool creation outperforms prior agentic pipelines, with the largest gains on complex and domain\-specific tasks\. We further show that dynamically generated tools are substantially more sophisticated than those produced by existing visual tool\-creation systems, covering more diagnostic categories per function and strongly preferring visual outputs that the VLM critic can reason over directly\.

## 1Introduction

Large Language Model \(LLM\) and Vision\-Language Model \(VLM\) agents have demonstrated remarkable promise in scientific discovery, with systems capable of forming hypotheses, designing experiments, and synthesizing results\. However, failures at different steps impede these models’ abilities to truly accelerate science\. This work focuses on the task offitting quantitative models to data\. This step is complex, requiring priors about the domain in question and careful statistical analysis, while being less studied than other quantitative analysis tasks like ML engineering and data science\[[13](https://arxiv.org/html/2606.00384#bib.bib46),[10](https://arxiv.org/html/2606.00384#bib.bib47)\]\. Addressing this gap is particularly important in scientific disciplines like astronomy where models struggle to effectively use domain conventions\[[23](https://arxiv.org/html/2606.00384#bib.bib55)\]\.

We introduceVesta111VESTA is the second most massive body in the main asteroid belt\[[33](https://arxiv.org/html/2606.00384#bib.bib50)\]\.\(VisualExploration withStatisticalToolAgents\), a framework that proposes, refines and evaluates statistical models by analyzingvisual representationsof data\.Vestaoperates over models defined in PyMC\[[1](https://arxiv.org/html/2606.00384#bib.bib53)\], a Bayesian probabilistic programming framework that enables model definition and fitting to data using an efficient Markov Chain Monte Carlo \(MCMC\) method\. Model\-building with PyMC code is compositional, allowing for complex distributions to be constructed by combining simpler components\.Vestainstantiates a loop of proposing models, selecting the best one, then iterating, broadly similar to existing autoresearch workflows\[[40](https://arxiv.org/html/2606.00384#bib.bib49),[8](https://arxiv.org/html/2606.00384#bib.bib48)\]\. However, unlike prior work,Vestautilizes multimodal feedback by augmenting this basic loop with a new type ofvisual tools\. Visual tools allowVestato inspect model fits using a VLM, giving it leverage in a modality distinct from past work on tool use\.

Figure[1](https://arxiv.org/html/2606.00384#S1.F1)shows how a model can use a tool to compute residuals, produce a visualization, and then inspect that visualization to identify areas where the predicted model does not fit the data well and propose a new statistical model that better fits the data\. For the diversity of statistical model fitting tasks we explore here, it is infeasible to enumerate all necessary tools up front;*dynamic*tool instantiation is required to cover the full breadth of cases\.

![Refer to caption](https://arxiv.org/html/2606.00384v1/x1.png)Figure 1:Overview ofVesta\. By effectively using and creating tools,Vestaproduces a probabilistic PyMC program that models the input data\.To support advancement on challenging model fitting problems where current approaches fail, we introduce a new benchmark we callDawn\(Dataset forAutomatedWorkflows andNumerical Modeling\)\.222DAWN was NASA’s first truly interplanetary spaceship that explored VESTA\[[32](https://arxiv.org/html/2606.00384#bib.bib51)\]\.Dawnproblems cover two domains central to both general data science and scientific research: fitting probability distributions to discrete data and constructing models for time\-series analysis\. In contrast to existing benchmarks which often contain a small number of distributions with a limited number of points\[[15](https://arxiv.org/html/2606.00384#bib.bib12),[29](https://arxiv.org/html/2606.00384#bib.bib54)\],Dawnconsists of synthetic problems spanning two levels of difficulty \(EasyandHard\), as well as two fundamental problems from astronomy \(Astro\): fitting stellar initial mass functions \(IMFs\), and detecting and characterizing chirp signals embedded in gravitational\-wave data\.

UsingDawn, we benchmarkVestaagainst prior frameworks specialized for Bayesian model and dynamic VLM\-based tool\-generation systems\. While existing approaches perform competitively on ourEasysplits, performance degrades substantially as data complexity increases and tasks shift toward domain\-specific scientific phenomena\. This degradation reveals a systematic gap in current agents’ capacity for iterative, evidence\-driven model refinement\. We show that this is effectively addressed byVesta’s visual tools, with which the agent is capable of exploring multiple hypotheses, visualizing them, and generating new hypotheses based on its tool’s visual analyses\. Furthermore, we find thatVestanot only successfully recovers the expert tools at a higher rate than existing baselines, but also creates sophisticated tools that often composes several expert tools in order to test multiple hypotheses at once\. While using dynamically generated visual tools does not outperform expert\-written tools, likely due to limitations in agents’ ability to process complex diagrams,Vestashows great promise for creating and utilizing statistical tools\.

## 2Related Work

##### Automated Model Fitting

Past work evaluated the capabilities of LLMs and VLMs to automate statistical analysis\[[42](https://arxiv.org/html/2606.00384#bib.bib9),[22](https://arxiv.org/html/2606.00384#bib.bib11)\], time\-series modeling\[[52](https://arxiv.org/html/2606.00384#bib.bib5),[47](https://arxiv.org/html/2606.00384#bib.bib4),[16](https://arxiv.org/html/2606.00384#bib.bib3),[56](https://arxiv.org/html/2606.00384#bib.bib6),[34](https://arxiv.org/html/2606.00384#bib.bib8)\]and Bayesian model proposals\[[37](https://arxiv.org/html/2606.00384#bib.bib19),[28](https://arxiv.org/html/2606.00384#bib.bib15)\]\.Jiet al\.\[[22](https://arxiv.org/html/2606.00384#bib.bib11)\]andSunet al\.\[[42](https://arxiv.org/html/2606.00384#bib.bib9)\]survey how LLM\-based agents can automate statistical workflows from data cleaning and model selection\. However,Zhuet al\.\[[59](https://arxiv.org/html/2606.00384#bib.bib10)\]find that even state\-of\-the\-art LLMs achieve only modest accuracy on benchmarks for selecting appropriate statistical methods, revealing significant limitations in LLMs’ ability to reason over complex statistical tasks\. While much of the focus in the literature has used LLMs / VLMs to predict a continuation of time series or detect anomalies, our work seeks to representunderlying modelthat generates the observed data using a probabilistic Bayesian model \(PyMC\)\. Previous work in generating PyMC code to model data has been successful in scientific settings, but this work often relies on human\-written kernels for Gaussian process, limiting the use cases of such methods\[[5](https://arxiv.org/html/2606.00384#bib.bib18),[6](https://arxiv.org/html/2606.00384#bib.bib16),[12](https://arxiv.org/html/2606.00384#bib.bib17)\]\. Recent works addresses these key limitation by using LLMs to automatically perform “Box’s Loop” of iteratively proposing new and critiquing existing models written in PyMC code until an adequate fit is achieved\[[27](https://arxiv.org/html/2606.00384#bib.bib13)\]\. Despite some promising results using LLMs to propose probabilistic models,Gandhiet al\.\[[15](https://arxiv.org/html/2606.00384#bib.bib12)\]find LLMs struggle on more challenging datasets sampled from real\-world domains and demonstrate that the iterative refinement of PyMC models does not consistently provide better predictions over one\-shot LLMs\. Further, while existing benchmarks sample data across different domains, they cover only a narrow range of distributions, limiting their difficulty and scope\. OurVestaautomates Box’s propose\-critique\-refine loop using VLMs, enabling tool use and creation for data analysis, while handling significantly more complex distributions\.

##### Multimodal Tool Use & Creation

Several works have demonstrated that VLMs are proficient tool\-users and that the successful utilization of tools can improve fine\-grained reasoning over visual\-question answering tasks\[[20](https://arxiv.org/html/2606.00384#bib.bib22),[57](https://arxiv.org/html/2606.00384#bib.bib23),[49](https://arxiv.org/html/2606.00384#bib.bib24),[51](https://arxiv.org/html/2606.00384#bib.bib28)\]\. Incorporating VLM tool\-use in agentic systems has produced even larger performance gains\[[53](https://arxiv.org/html/2606.00384#bib.bib29),[18](https://arxiv.org/html/2606.00384#bib.bib30),[50](https://arxiv.org/html/2606.00384#bib.bib31),[19](https://arxiv.org/html/2606.00384#bib.bib32)\]\. In particular,Zhanget al\.\[[54](https://arxiv.org/html/2606.00384#bib.bib21)\]propose a multi\-agent VLM framework in which an orchestrator agent coordinates specialized sub\-agents and vision expert tools to improve fine\-grained visual perception\. Recent work has tested the capabilities of VLMs; for example, in PyVision\[[55](https://arxiv.org/html/2606.00384#bib.bib20)\], VLMs are used in an agentic flow there they dynamically create and execute tools until an answer is reached\. Previous work has limited VLM tool use \(i\.e\., executing or generating code over visual inputs\) to image enhancements, annotations, cropping, and the calculation of basic statistics\[[17](https://arxiv.org/html/2606.00384#bib.bib25),[44](https://arxiv.org/html/2606.00384#bib.bib26),[58](https://arxiv.org/html/2606.00384#bib.bib27)\]\. While some work has applied multimodal tool creation to model fitting\[[24](https://arxiv.org/html/2606.00384#bib.bib14)\], such tools are limited to simple statistics and basic plotting, lack a dedicated tool\-creation step, and are discarded after each use rather than persisted for reuse across iterations\. In contrast,Vestauses an LLM tool agent to dynamically create visualizations that test statistical hypotheses, which are then passed to a VLM critic as evidence for refining them\.

##### Autonomous Tool Expansion

While both LLMs and VLMs have proven to use tools effectively, curating the exact set of tools for a problem often requires manually curating a tool kit\. Initial work demonstrates that LLMs are capable of coding effective and reliable tools\[[36](https://arxiv.org/html/2606.00384#bib.bib34),[7](https://arxiv.org/html/2606.00384#bib.bib33)\], although these tools are typically created separately from model inference\[[2](https://arxiv.org/html/2606.00384#bib.bib36),[48](https://arxiv.org/html/2606.00384#bib.bib35)\]\. Systems such as Voyager\[[46](https://arxiv.org/html/2606.00384#bib.bib37)\]extend this static paradigm by allowing for tool libraries to adapt to task needs and evolve across examples\. Many of the systems for autonomous tool expansion tend to focus on software engineering tasks\[[11](https://arxiv.org/html/2606.00384#bib.bib39),[43](https://arxiv.org/html/2606.00384#bib.bib40),[26](https://arxiv.org/html/2606.00384#bib.bib45)\]or narrowly scoped API invocation benchmarks\[[21](https://arxiv.org/html/2606.00384#bib.bib38),[41](https://arxiv.org/html/2606.00384#bib.bib43)\], with less focus on problems that require, judgment\-driven statistical analysis\. While EvoSkill\[[2](https://arxiv.org/html/2606.00384#bib.bib36)\]and TTE\[[30](https://arxiv.org/html/2606.00384#bib.bib41)\]create some tools related to calculating statistics, they tend to contain more simplistic functions such as calculating averages, areas or retrieving information from tables\.

## 3Vesta: Visual Exploration with Statistical Tools

We assume as input a datasetDDgenerated by a modelM​\(θ\)M\(\\theta\)\.MMis a data generating process \(e\.g\., a directed graphical model\) with parametersθ\\theta\. In this work, we consider models that can be specified in PyMC code\. Given a proposed modelM~​\(θ~\)\\tilde\{M\}\(\\tilde\{\\theta\}\), we can evaluate its fit toDDwith metrics denoted asR​\(D,M~​\(θ~\)\)R\(D,\\tilde\{M\}\(\\tilde\{\\theta\}\)\)\. Central to our approach arevisual tools,ε​\(D,M,θ\)\\varepsilon\(D,M,\\theta\), implemented as executable Python functions that accept the dataDDand a current model\(M,θ\)\(M,\\theta\)and return visual or statistical diagnostics used to inform subsequent modeling decisions\.

### 3\.1Vesta

Algorithm[1](https://arxiv.org/html/2606.00384#alg1)illustrates the workflow forVesta\. We start with an initial prediction by estimating a probabilistic model from an initial plot of the data, and then refine this prediction overNNsteps\. At each step, we proposeppalternative candidate models, then select model with the best metric,RR\. Our approach proceeds in four steps per iteration: propose, generate tools, critique model fit, and summarize steps of the run\.

Propose\.On the first iteration,Vestatakes an initial visualization of the data, produces a brief visual description and is instructed to createppdiverse model proposals\. To ensure faithfulness to the proposed model parameters, PyMC code is generated in a dedicated step that translates parameter specifications into executable PyMC programs\. Without this step, we find models frequently deviate from their specified parameterizations\.

Algorithm 1Visual Exploration Agents1:Data

DD, iterations

NN, proposals

pp, metric

RR, registry

ℰ\\mathcal\{E\}\(initial state:generate\_new\_toolonly\)

2:

Mbest,θbestM\_\{\\text\{best\}\},\\theta\_\{\\text\{best\}\}
3:

Mbest,θbest←Propose​\(Plot​\(D\)\)M\_\{\\text\{best\}\},\\theta\_\{\\text\{best\}\}\\leftarrow\\textsc\{Propose\}\(\\text\{Plot\}\(D\)\)⊳\\trianglerightPlot data; fitjjmodels by and select best model under metric,RR

4:

s0←Summarize​\(Mbest,θbest\)s\_\{0\}\\leftarrow\\textsc\{Summarize\}\(M\_\{\\text\{best\}\},\\theta\_\{\\text\{best\}\}\)
5:for

i=1,…,Ni=1,\\ldots,Ndo

6:

ε∗←ToolManager​\(si−1,D,ℰ\)\\varepsilon^\{\*\}\\leftarrow\\textsc\{ToolManager\}\(s\_\{i\-1\},D,\\mathcal\{E\}\)⊳\\trianglerightSelect/create tool; updateℰ\\mathcal\{E\}

7:

Ti←ε∗​\(D,Mi−1,θi−1\)T\_\{i\}\\leftarrow\\varepsilon^\{\*\}\(D,M\_\{i\-1\},\\theta\_\{i\-1\}\)⊳\\trianglerightExecute tool

8:

Mi,θi←Critique​\(Ti,s0:i−1,D\)M\_\{i\},\\theta\_\{i\}\\leftarrow\\textsc\{Critique\}\(T\_\{i\},s\_\{0:i\-1\},D\)⊳\\trianglerightProposeppmodels; return best ranked byRR

9:

Mbest,θbest←arg⁡minθ∈\{θi,θbest\}⁡R​\(ℳ​\(θ\)\)M\_\{\\text\{best\}\},\\theta\_\{\\text\{best\}\}\\leftarrow\\arg\\min\_\{\\theta\\in\\\{\\theta\_\{i\},\\theta\_\{\\text\{best\}\}\\\}\}R\(\\mathcal\{M\}\(\\theta\)\)
10:

si←Summarize​\(Mi,θi\)s\_\{i\}\\leftarrow\\textsc\{Summarize\}\(M\_\{i\},\\theta\_\{i\}\)
11:endfor

12:return

Mbest,θbestM\_\{\\text\{best\}\},\\theta\_\{\\text\{best\}\}

Generate Tools\.The central contribution ofVestais its ability to explore data through the composition and creation of tools\. At each iteration,Vestareceives a summary of the previous runs and model parameters then decides to either selects an existing tool from a growing registryℰ\\mathcal\{E\}, or dynamically codes a new diagnostic tool in Python when no existing tool is adequate for the desired analysis\. Newly created tools are saved to the registry, making them reusable across future iterations for the current run\. At each run, the toolkit is newly initialized\. This allowsVestato construct increasingly sophisticated, hypothesis\-driven visualizations\. We find in Section[6\.2](https://arxiv.org/html/2606.00384#S6.SS2)that tools created byVestahave high coverage of expert written tools and create sophisticated visualizations to test multiple statistical hypotheses in a single tool call\.

Critique\.The output of the selected or newly created tool, along with the history of the best models and descriptions from previous runs is passed to a VLM critic that refines the current PyMC model by proposingpprevised models\. The revised model that optimizesRRis passed to the next iteration\. In our prompt, we instruct the VLM critic to terminate the iterative refinement loop early if it judges the current model fit to be satisfactory\. This is implemented by instructing the critic to include an explicitCOMPLETEtoken in its description, which causesVestato return the best PyMC code\.

Summarize\.At the end of each iteration, the visual description, chosen model and parameters, the metricRR, and tool output are compressed into a structured summary\. This summary is injected into the next iteration in place of raw intermediate outputs, enabling the system to reason over the full refinement trajectory without incurring prohibitive context growth\.

### 3\.2Visual Tools inVesta

Vestacan work out\-of\-the\-box without tools, a setting we evaluate\. However, we found tools to be critical in solving challenging model fitting tasks\. In particular, we allowVestato propose new tools tailored to their model predictions\. At each tool generation step,Vestareceives as input the summarized contents of the previous iteration \(including the VLM critique of the model fit, the best PyMC model withR​\(D,Mbest​\(θ\)\)R\(D,M\_\{\\text\{best\}\}\(\\theta\)\), and a history of previously executed tools, along with access to a tool registry\. We initialize the toolkit with a single function that gives the option for creating a new tool\. New tools emerge through the process of examining data, and allow models to verify, and update an existing hypothesis \(e\.g\., a proposed statistical distribution\)\. As we progress through the task, we save previously created tools\. If a generated tool fails, we allow up to two retries passing the failure trace into the retry prompt\. After two retries if a tool fails to execute due to syntax or import errors, the system falls back to the “no toolkit” mode\.

In practice, we found thatVestais capable of generating complex tools that tests distinct hypotheses through visualization, which is then utilized in the Critique stage\. The tools always involve visual transformations to explore the data, which we show in detail in a case study in Section[6\.2](https://arxiv.org/html/2606.00384#S6.SS2)\(Figure[5](https://arxiv.org/html/2606.00384#S6.F5)\)\.Vestacan also work with externally provided visual tools\. In Section[5](https://arxiv.org/html/2606.00384#S5)we compile a list ofexpert\-writtentools which we evaluate and analyze in Section[6](https://arxiv.org/html/2606.00384#S6)\.

## 4TheDawnBenchmark

Distribution fitting and time series modeling are two key data science modeling tasks that appear consistently across scientific disciplines\. We select these domains because they allow us to benchmark AI agents across easily recognizable models, complex mixtures of common forms, and challenging real\-world tasks where domain\-specific knowledge, such as established conventions and governing equations, can reduce the complexity of modeling real\-world phenomena\.Dawnproblems are synthetically generated, allowing us to robustly measure model fit\. Figure[2](https://arxiv.org/html/2606.00384#S4.F2)shows examples for each domain and split \(Easy,Hard, andAstro\)\.

![Refer to caption](https://arxiv.org/html/2606.00384v1/images/data_ex.png)Figure 2:Sample inputs from both domains and all dataset splits inDawn\.Easysplits contain easily recognizable forms\.Hardtasks contain mixtures of distinct forms, andAstrotasks reflect real\-world astronomy challenges that require additional analysis beyond simple visualization to solve\.### 4\.1Distribution Fitting Tasks

TheEasytasks in distribution fitting consists of identifying the family and associated parameters for a unimodal distribution\. The distribution families we consider in this paper are:gaussian,lognormal,student\-t,exponential,uniform,weibull,laplace,cauchyandpareto\. Both the priors and exact number of points for each distribution are randomly sampled\. TheHardtask composes exactly two of these distributions\. Mixtures present in theHardsplit model phenomenon present other real\-wold domains such as the gaussian\-laplace mixtures common in speech signal processing\[[3](https://arxiv.org/html/2606.00384#bib.bib62)\]and lognormal\-exponential mixtures prevalent in survival\-analysis\[[31](https://arxiv.org/html/2606.00384#bib.bib63)\]\.

For ourAstrosplit, we use initial mass functions \(IMF\) as the distribution family, which describes the distribution of stellar masses at birth and are often modeled with canonical forms\[[35](https://arxiv.org/html/2606.00384#bib.bib56)\]\. Accurate modeling of the IMF allows astronomers to infer the conditions of past star formation and predict its downstream effects on galactic evolution and heavy element distribution\[[4](https://arxiv.org/html/2606.00384#bib.bib52)\]\.

We use five functional forms \(visualized in Figure[3](https://arxiv.org/html/2606.00384#S4.F3)\):\(1\)The simplest form, the Salpeter IMF\[[39](https://arxiv.org/html/2606.00384#bib.bib57)\], is a single power law\.\(2\)The Kroupa IMF\[[25](https://arxiv.org/html/2606.00384#bib.bib58)\]is a piecewise power law with sharp breaks at specific characteristic masses\.\(3\)The Chabrier IMF\[[9](https://arxiv.org/html/2606.00384#bib.bib59)\]is distinguished by a smooth log\-normal component at low masses transitioning to a power law at higher masses\. We additionally include two freeform variants that do not follow canonical forms making them more challenging to model:\(4\)one with breaks constrained near canonical masses \(freeform tight\) and\(5\)one with unconstrained breaks and slopes that can produce more exotic stellar IMFs \(freeform wide\)\.

![Refer to caption](https://arxiv.org/html/2606.00384v1/images/imf_ex.png)Figure 3:Dawn’sAstrodistribution fitting tasks\. Example of Initial Mass Functions projected into log\-log space\. Distributions become visually distinct only when projected into log\-log space\.##### Evaluation Metric\.

We use Jensen\-Shannon Divergence to measure similarity between our proposed model’s distribution and the reference model’s distribution\. This metric is symmetric and bounded in the interval\[0,1\]\[0,1\], with a score of 0 indicating identical distributions\.

### 4\.2Time Series Fitting Tasks

TheEasytask in time series fitting consists of modeling sequences characterized by simple linear trends \(increasing, decreasing, or flat\) and standard seasonal components, such as basic periodic sine waves, coupled with minimal noise\. These datasets are designed to evaluate whether models can correctly identify and combine standard components such as linear and periodic kernels with minimal ambiguity\. TheHardtask introduces more complex, non\-linear dynamics and non\-standard periodicities\. This split includes time series exhibiting random walks \(ARIMA process\), and S\-curve \(sigmoidal\) functions\. Further,Hardcontains time series from other domains, such as electrocardiogram \(ECG\) signals of heartbeat patterns and time series exhibiting seasonal periodicity with an increasing averages common in sales forecasting\[[14](https://arxiv.org/html/2606.00384#bib.bib64)\]\. Capturing the behavior of these series requires moving beyond standard linear or single\-period functions to employ multi\-periodic modeling and Gaussian process kernels, such as Matérn and Radial Basis Function \(RBF\)\. This split tests the model’s ability to handle richer structures and moderate levels of noise and non\-linearity\.

For ourAstrosplit, we model time series inspired by gravitational wave chirps\. “Chirps” are an astronomical phenomenon produced by merging binary star systems that introduces severe non\-stationarity into the signal\. The standard chirp variant features a continuously changing frequency over time, where initial high\-wavelength, low\-frequency waves progressively compress into higher frequencies as the celestial bodies converge\. We additionally include a highly complex variant where this dynamic, shifting frequency is coupled with an amplitude that decays over time\.

##### Evaluation Metric\.

In accordance with prior works\[[27](https://arxiv.org/html/2606.00384#bib.bib13)\], we use the Expected Log Predictive Density under Leave One Out \(ELPD\-LOO\) metric for time series\. Formally, assume observationsy1,…,yny\_\{1\},\\dots,y\_\{n\}are independent given parametersθ\\theta\. Letp​\(θ\)p\(\\theta\)be a prior distribution on the model parameters that forms a posterior distributionp​\(θ∣y\)p\(\\theta\\mid y\)and posterior predictivep​\(y~∣y\)=∫p​\(y~∣θ\)​p​\(θ∣y\)​𝑑θp\(\\tilde\{y\}\\mid y\)=\\int p\(\\tilde\{y\}\\mid\\theta\)p\(\\theta\\mid y\)\\,d\\theta\. We measure predictive accuracy using the expected log pointwise predictive density:

ELPD−LOO=∑i=1n∫pt​\(y~i\)​log⁡p​\(y~i∣y\)​𝑑y~i\\mathrm\{ELPD\-LOO\}=\\sum\_\{i=1\}^\{n\}\\int p\_\{t\}\(\\tilde\{y\}\_\{i\}\)\\,\\log p\(\\tilde\{y\}\_\{i\}\\mid y\)\\,d\\tilde\{y\}\_\{i\}\(1\)wherept​\(y~\)p\_\{t\}\(\\tilde\{y\}\)denotes the true data\-generating distribution that we use LOO to approximate\. ELPD measures predictive accuracy on a log\-probability scale, with higher values indicating better fit\. While unbounded below, its absolute value is less important than differences between models, which quantify relative predictive performance\. Due to the computational cost of running ELPD\-LOO over the full time series, we apply a sliding\-window sampling procedure to create a representative subsample \(25%\) of the signal\.

### 4\.3Dataset Setup and Statistics

LetD¯=\{D1,…,Dn\}\\bar\{D\}=\\\{D\_\{1\},\\ldots,D\_\{n\}\\\}be a collection of datasets where eachDiD\_\{i\}represents a set ofn∈\[600,1500\]n\\in\[600,1500\]points randomly sampled from a ground\-truth probability distribution or time\-series signalDTrueiD\_\{\\text\{True\}\_\{i\}\}, whose parameters are drawn uniformly from fixed ranges \(see Appendix[B](https://arxiv.org/html/2606.00384#A2)\)\. For each domain and difficulty split we generate 50 \(Easy/Astro\) or 100 \(Hard\) such datasets, evaluated against the known ground\-truth model\. In total, our benchmark consists of 200 discrete probability distributions and 200 time series\.

## 5Experimental Setup

### 5\.1VestaImplementation Details

We use the following models as backbones for ourVesta: GPT\-5\.4\-mini, Claude Sonnet 4\.6, and Kimi K2\.5\. Note that the backbone model forVestamust support interleaved text and image input\. AlthoughVestacan be implemented with different models at each step, we use the same one for simplicity\. For GPT\-5\.4\-mini the reasoning effort is set to “low” and for Claude\-Sonnet\-4\.6 we select the minimal number of reasoning tokens \(1024\) due to the iterative nature of our task and exploding token lengths\.

For all experiments run withVestawe set the number of max iterationsN=5N=5, number of proposals per stepp=3p=3, and we selectRRto be the Akaike Information Criterion \(AIC\) as it balances model complexity with data fit and does not depend on having access to a ground truth distribution\. For the initial plot, we use a histogram for distributions and simple line plot for time series \(See Figure[2](https://arxiv.org/html/2606.00384#S4.F2)\)\.

##### Expert Visual Tools

In addition to dynamically generating tools, we can also equipVestawith tools designed by an expert statistician as an*oracle*experiment\. The expert toolkit remain unaltered during the course of exploration and consists of general\-purpose, functions that target key questions expert statisticians use when fitting complex mixture distributions or noisy time\-series\. For distribution fitting, the static toolkit includes tools for visualizing QQ plots, tail transform diagnostics, probability plots, and a running a Gaussian\-Mixture\-Model based segmentation with moment estimation, enabling agents to systematically diagnose location, scale, tail, and mixture structure miscalculations\. For time series, the expert toolkit includes tools for comparing model fit against actuals, diagnosing residuals, and testing temporal independence via autocorrelation plots and Ljung\-Box statistics, together enabling agents to detect underfitting, missed periodicity, and unmodeled dynamics\. A complete description of tools is available in Appendix[E](https://arxiv.org/html/2606.00384#A5)\.

### 5\.2Baselines

We detail exact prompts of all baselines in Appendix[H](https://arxiv.org/html/2606.00384#A8)\. For fair comparison withVesta, we set the maximum number of iterationsN=5N=5for all baselines\.

##### BoxLM\.

BoxLM\[[27](https://arxiv.org/html/2606.00384#bib.bib13)\]is inspired by the principle of Box’s Loop where a language model proposes statistical models using PyMC code and then critiques the probabilistic program overNNrounds of iteration\. The LLM proposesMMmodels per round and scores each of them using ELPD\-LOO, selecting the topppas candidates for the next round\. The critic LM then synthesizes the fitted models, their scores, and posterior predictive summary statistics into a natural language hypothesis in order to update the model proposals in the next round of feedback\. We minimally adapt Box LM prompts to better reflect our task goals\.

##### PyVision\.

PyVision\[[55](https://arxiv.org/html/2606.00384#bib.bib20)\]is a multi\-turn framework designed to enhance multimodal reasoning by allowing VLMs to generate and execute Python code during inference\. A VLM receives an input, generates Python code, and executes it in an isolated runtime, with the resulting output fed back into the VLM’s context to refine reasoning over multiple turns\. The iteration continues until the VLM decides to provide a final answer\. PyVision is originally designed for Visual Question Answering, where it generates and executes tools to enhance images by adjusting contrast, cropping irrelevant regions, and computing simple statistics such as histograms of pixel intensities\. We adapt PyVision to propose and fit PyMC models by providing the same intent and PyMC coding guidelines that Visual Exploration Agent receives, keeping the multi\-turn execution loop intact\.

## 6Results and Analysis onDawn

![Refer to caption](https://arxiv.org/html/2606.00384v1/images/js_elpd_bar.png)Figure 4:\[Top\]Average Jensen\-Shannon divergence \(↓\\downarrowbetter\) between the ground\-truth distribution and the probability density function of the proposedPyMCmodel on the Distribution Fitting task ofDawn\.\[Bottom\]Average ELPD\-LOO \(↑\\uparrowbetter\) for the Time Series Modeling task ofDawn, computed via leave\-one\-out cross\-validation\. Error bars denote±\\pm1 standard error of the mean\.### 6\.1Results

##### Distribution Fitting

Figure[4](https://arxiv.org/html/2606.00384#S6.F4)reports the Jensen\-Shannon Divergence between the ground truth distribution and the probability density function induced by agent\-written PyMC code\.Vestaequipped with tool use consistently outperforms both the PyVision and BoxLM baselines acrossDawn, with the largest gains on theHardandAstrosplits\. When provided with an expert toolkit,Vestaachieves the strongest overall performance, establishing an upper bound for what tool discovery pipelines can attain\. Crucially,Vestawith dynamic tool generation approaches this upper bound without access to any expert\-written tools, producing distributions that far more closely reflect ground truth structure than the no\-tool baseline\. This gap demonstrates thatVestais an effective tool discovery method: autonomously generated tools are sufficient to capture the modeling inferences that expert tools provide\.

##### Time Series

Similar to the results for Distribution Fitting,Vesta, with both dynamic toolkit and expert toolkit consistently beats both PyVision and BoxLM baseline across theDawnTime Series data\.Vestaequipped with the expert toolkit has the strongest performance, followed closely byVestawith the dynamic toolkit\. We find that current baselines especially struggle to effectively model the time series inDawn\.

### 6\.2Analysis of Generated Tools

![Refer to caption](https://arxiv.org/html/2606.00384v1/images/vesta_ex.png)Figure 5:Example of the output from aVestagenerated tool\. This tool composes multiple functions to analyze a heavy\-tailed distribution\. This multi\-panel visualization output is fed back intoVestato generate better hypotheses\. Panel titles are enlarged for clarity and panel numbers are added manually\.##### Case Study

Figure[5](https://arxiv.org/html/2606.00384#S6.F5)shows an example of aVestawritten tool designed to determine the family of a given right\-tailed distribution\. The tool generated a single visualization with 6 panels; the original histogram of the input distribution is inPanel 4\. Each of the remaining subplots are designed to test particular hypothesis about this skewed distribution\.Panel 1shows a Q\-Q plot testing whether the current proposed exponential distribution based PyMC model is an appropriate for this data, shown by whether points fall along the reference line\.Panel 2plots a histogram of the raw data on a log\-scaledxx\-axis, which compresses the right tail and spreads out the lower values, making it easier to assess the shape of the distribution and identify whether a heavy tail is present\.Panel 3plots a histogram of the log\-transformed data, meaning each observation is logged before being binned\. If this histogram appears bell\-shaped and symmetric then a log\-normal family is appropriate for the original data\.Panel 5is a box plot andPanel 6plots the complementary cumulative distribution function \(CCDF\) of the data in log\-log space to see if the data follows a power\-law distribution\. This multi\-panel visualization is a very powerful ingredient forVesta; one limitation is that VLMs sometimes exhibit a perceptual “blindness”\[[38](https://arxiv.org/html/2606.00384#bib.bib60),[45](https://arxiv.org/html/2606.00384#bib.bib61)\]and fail to analyze these complex visualizations, a technical limitation that if addressed, will likely yield further performance gains\.

##### Generated vs\. Expert Tools

We analyze the tools generated byVestaand PyVision using Claude Sonnet 4\.6, and compare them to the expert\-written toolkit\. We find thatVestaand PyVision often independently “discover” tools present in the expert toolkit and compose them together to create a single, powerful tool for analysis\. The tool in Figure[5](https://arxiv.org/html/2606.00384#S6.F5)provides detailed analysis in a single tool that would require multiple steps\. In contrast, expert tools are atomic by design, often testing a maximum of two hypotheses\. For example, a single tool generated byVestacomposes an average of 3\.5 tools present in the expert toolkit, compared to 1\.9 tools for PyVision\. This finding indicates thatVestais capable of creating powerful statistical tools that have a high overlap with tools written by human experts, and often explore multiple hypotheses simultaneously\.

Table[1](https://arxiv.org/html/2606.00384#S6.T1)quantifies this analysis at scale\. We calculate the percent of times a generated tool contains a sub\-function present in the expert toolkit\. First, we find that both methods discover each function type in the tool kit at least once\.Vestarecovers a much larger percentage of expert tools than PyVision\.

Second, we find that generated tools depend on the unique specifications of the task\. For distribution fitting, theAstrosplit contains Initial Mass Functions \(IMFs\) which are exclusively power law or mixtures of power law distributions\. In order to visually distinguish these distributions, applying a log or log\-log transform is required\. Accordingly, forVestathe percent of generated tools using a “tail transform” function increases from 55\.4% on theEasysplit to 91\.8% on theAstro\.

Table 1:Expert tool coverage \(%\) across difficulty levels forVestaand PyVision using Claude Sonnet 4\.6\. Each value is the percentage of generated tools that implement one of the expert tools\.
##### VestaCreated Tools Beyond the Expert Toolkit

Not only doesVestadiscover expert\-written tools, but it creates tools that further strengthen analysis\.

Table 2:For the distribution fitting task, the percentage ofVesta\-generated tools that invoke each function as a subroutine\. Only functions not present in the expert toolkit are shown\.Table[2](https://arxiv.org/html/2606.00384#S6.T2)shows the percentage ofVesta\-generated tools that contain a function not present in the expert toolkit\. We find that nearly all generated tools contain some form of “diagnostic check” where the tool fits one more candidate distributions against the input data to assess which distribution best models\. This is often accompanied by a formal ranking by calculation different “information criterion” metrics such as AIC or BIC\. We find the next most common functions \(mean\-excess and Hill estimator\) are predominantly called on theAstrosplit and are designed to diagnose heavy\-tailed behavior\. The mean\-excess plot exploits linearity of𝔼​\[X−u​∣X\>​u\]\\mathbb\{E\}\[X\-u\\mid X\>u\]in the thresholduuto identify a Generalized\-Pareto\-like tail and locate where it begins, while the Hill estimator provides a direct estimate of the tail indexα\\alphafrom the largest order statistics\. Their focus onAstrosuggests thatVestacorrectly identifies Initial Mass Functions as tail\-dominated distributions and tailors its tools to target these features\. The remaining tools \(Box–Cox, D’Agostino’s normality test, and Shapiro–Wilk\) are all related to Gaussian distributions, either testing for normality or transforming data to be normally distributed\. A detailed description of each tool, as well as an analysis of novel time\-series tools present in Appendix[F](https://arxiv.org/html/2606.00384#A6)\.

##### Tool Persistence and Reuse\.

A key architectural advantage of VESTA’s dynamic\-toolkit is the persistence of its generated functions\. Within a single problem instance, 79\.5% of VESTA’s dynamically generated tools are re\-invoked across subsequent analysis steps \(averaging 1\.87 calls per tool; Table[3](https://arxiv.org/html/2606.00384#S6.T3)\)\. This capability is most useful for tools that evaluate model fit, allowing the agent to iteratively reassess model fit using a consistent baseline\. Consider the case study shown in Figure[5](https://arxiv.org/html/2606.00384#S6.F5)\. AsVestaiteratively refines its distributional hypothesis \(from exponential to log\-normal\), the tool can be re\-invoked with an updated PyMC model to regenerate Panel 1’s Q\-Q plot against the new candidate distribution without incurring additional tool\-writing cost\. Further, the remaining diagnostic panels act as a consistent reference frame, meaningVestadoes not have to regenerate tools to re\-examine the underlying data when evaluating a new distributional hypothesis\. PyVision treats generated code as single\-use execution cells, entirely lacking a mechanism for revisitation\.

Table 3:VESTA tool reuse statistics within single problem instances\.

## 7Conclusion

We introduceVesta, a pipeline that uses VLMs to propose statistical models and test hypotheses through dynamic tool use, producing models that more accurately reflect underlying data structure\. To support this evaluation, we introduceDawn, a benchmark of distribution fitting and time series modeling tasks at varying difficulty, including a challengingAstrosplit drawn from real\-world astronomy problems\. Our evaluation shows thatVestaindependently recovers expert tools and composes them into more sophisticated diagnostics that test multiple hypotheses simultaneously\. WhileVestawith expert tools achieves the strongest overall performance,Vesta\-Dynamic outperforms all existing baselines onDawn, demonstrating thatVestacan propose high\-quality statistical models and perform autonomous statistical tool discovery\.

## Acknowledgments

This work was supported by the NSF under Cooperative Agreement 2421782 and the Simons Foundation grant MPS\-AI\-00010515 awarded to the NSF\-Simons AI Institute for Cosmic Origins — CosmicAI,[https://www\.cosmicai\.org/](https://www.cosmicai.org/)\. It was also partially supported by NSF grants IIS\-2145280 and IIS\-2433071\. We would also like to thank Stéphanie Juneau from CosmicAI / NOIRLab for her help and feedback on this work\.

## References

- \[1\]O\. Abril\-Pla, V\. Andreani, C\. Carroll, L\. Dong, C\. J\. Fonnesbeck, M\. Kochurov, R\. Kumar, J\. Lao, C\. C\. Luhmann, O\. A\. Martin,et al\.\(2023\)PyMC: a modern, and comprehensive probabilistic programming framework in python\.PeerJ Computer Science9,pp\. e1516\.Cited by:[§1](https://arxiv.org/html/2606.00384#S1.p2.1)\.
- \[2\]\(2026\)EvoSkill: automated skill discovery for multi\-agent systems\.External Links:2603\.02766,[Link](https://arxiv.org/abs/2603.02766)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px3.p1.1)\.
- \[3\]A\. Aroudi, H\. Veisi, H\. Sameti, and Z\. Mafakheri\(2015\-12\)Speech signal modeling using multivariate distributions\.EURASIP Journal on Audio Speech and Music Processing2015,pp\. 1–14\.External Links:[Document](https://dx.doi.org/10.1186/s13636-015-0078-1)Cited by:[§4\.1](https://arxiv.org/html/2606.00384#S4.SS1.p1.3)\.
- \[4\]N\. Bastian, K\. R\. Covey, and M\. R\. Meyer\(2010\)A universal stellar initial mass function? a critical look at variations\.Annual Review of Astronomy and Astrophysics48\(Volume 48, 2010\),pp\. 339–389\.External Links:[Document](https://dx.doi.org/https%3A//doi.org/10.1146/annurev-astro-082708-101642),[Link](https://www.annualreviews.org/content/journals/10.1146/annurev-astro-082708-101642),ISSN 1545\-4282Cited by:[§4\.1](https://arxiv.org/html/2606.00384#S4.SS1.p2.1)\.
- \[5\]J\. Bongard and H\. Lipson\(2007\-06\)Automated reverse engineering of nonlinear dynamical systems\.Proceedings of the National Academy of Sciences of the United States of America104\(24\),pp\. 9943–9948\.External Links:[Document](https://dx.doi.org/10.1073/pnas.0609476104)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px1.p1.1)\.
- \[6\]J\. Brence, L\. Todorovski, and S\. Dzeroski\(2020\)Probabilistic grammars for equation discovery\.CoRRabs/2012\.00428\.External Links:[Link](https://arxiv.org/abs/2012.00428),2012\.00428Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px1.p1.1)\.
- \[7\]T\. Cai, X\. Wang, T\. Ma, X\. Chen, and D\. Zhou\(2024\)Large language models as tool makers\.External Links:2305\.17126,[Link](https://arxiv.org/abs/2305.17126)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px3.p1.1)\.
- \[8\]M\. Cemri, S\. Agrawal, A\. Gupta, S\. Liu, A\. Cheng, Q\. Mang, A\. Naren, L\. E\. Erdogan, K\. Sen, M\. Zaharia, A\. Dimakis, and I\. Stoica\(2026\)AdaEvolve: adaptive llm driven zeroth\-order optimization\.External Links:2602\.20133,[Link](https://arxiv.org/abs/2602.20133)Cited by:[§1](https://arxiv.org/html/2606.00384#S1.p2.1)\.
- \[9\]G\. Chabrier\(2003\-07\)Galactic stellar and substellar initial mass function\.Publications of the Astronomical Society of the Pacific115\(809\),pp\. 763–795\.External Links:ISSN 1538\-3873,[Link](http://dx.doi.org/10.1086/376392),[Document](https://dx.doi.org/10.1086/376392)Cited by:[§4\.1](https://arxiv.org/html/2606.00384#S4.SS1.p3.1)\.
- \[10\]J\. S\. Chan, N\. Chowdhury, O\. Jaffe, J\. Aung, D\. Sherburn, E\. Mays, G\. Starace, K\. Liu, L\. Maksin, T\. Patwardhan, L\. Weng, and A\. Mądry\(2025\)MLE\-bench: evaluating machine learning agents on machine learning engineering\.External Links:2410\.07095,[Link](https://arxiv.org/abs/2410.07095)Cited by:[§1](https://arxiv.org/html/2606.00384#S1.p1.1)\.
- \[11\]G\. Deng, Z\. Chen, Z\. Yu, H\. Fan, Y\. Liu, Y\. Yang, D\. Parikh, R\. Kannan, L\. Cong, M\. Wang, Q\. Zhang, V\. Prasanna, X\. Tang, and X\. Wang\(2026\)EvoClaw: evaluating ai agents on continuous software evolution\.External Links:2603\.13428,[Link](https://arxiv.org/abs/2603.13428)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px3.p1.1)\.
- \[12\]D\. Duvenaud, J\. R\. Lloyd, R\. Grosse, J\. B\. Tenenbaum, and Z\. Ghahramani\(2013\)Structure discovery in nonparametric regression through compositional kernel search\.External Links:1302\.4922,[Link](https://arxiv.org/abs/1302.4922)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px1.p1.1)\.
- \[13\]A\. Egg, M\. I\. Goyanes, F\. Kingma, A\. Mora, L\. von Werra, and T\. Wolf\(2025\)DABstep: data agent benchmark for multi\-step reasoning\.External Links:2506\.23719,[Link](https://arxiv.org/abs/2506.23719)Cited by:[§1](https://arxiv.org/html/2606.00384#S1.p1.1)\.
- \[14\]Y\. Ensafi, S\. Hassanzadeh Amin, G\. Zhang, and B\. Shah\(2022\)Time\-series forecasting of seasonal items sales using machine learning – a comparative analysis\.International Journal of Information Management Data Insights2\(1\),pp\. 100058\.External Links:ISSN 2667\-0968,[Document](https://dx.doi.org/10.1016/j.jjimei.2022.100058),[Link](https://www.sciencedirect.com/science/article/pii/S2667096822000027)Cited by:[§4\.2](https://arxiv.org/html/2606.00384#S4.SS2.p1.3)\.
- \[15\]K\. Gandhi, M\. Y\. Li, L\. Goodyear, A\. Bhatia, L\. Li, A\. Bhaskar, M\. Zaman, and N\. D\. Goodman\(2025\)BoxingGym: benchmarking progress in automated experimental design and model discovery\.External Links:2501\.01540,[Link](https://arxiv.org/abs/2501.01540)Cited by:[§1](https://arxiv.org/html/2606.00384#S1.p4.3),[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px1.p1.1)\.
- \[16\]N\. Gruver, M\. Finzi, S\. Qiu, and A\. G\. Wilson\(2024\)Large language models are zero\-shot time series forecasters\.External Links:2310\.07820,[Link](https://arxiv.org/abs/2310.07820)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px1.p1.1)\.
- \[17\]T\. Gupta and A\. Kembhavi\(2022\)Visual programming: compositional visual reasoning without training\.External Links:2211\.11559,[Link](https://arxiv.org/abs/2211.11559)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px2.p1.1)\.
- \[18\]J\. Hong, C\. Zhao, C\. Zhu, W\. Lu, G\. Xu, and X\. Yu\(2026\)DeepEyesV2: toward agentic multimodal model\.External Links:2511\.05271,[Link](https://arxiv.org/abs/2511.05271)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px2.p1.1)\.
- \[19\]X\. Hou, S\. Xu, M\. Biyani, M\. Li, J\. Liu, T\. C\. Hollon, and B\. Wang\(2026\)CodeV: code with images for faithful visual reasoning via tool\-aware policy optimization\.External Links:2511\.19661,[Link](https://arxiv.org/abs/2511.19661)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px2.p1.1)\.
- \[20\]Y\. Hu, W\. Shi, X\. Fu, D\. Roth, M\. Ostendorf, L\. Zettlemoyer, N\. A\. Smith, and R\. Krishna\(2024\)Visual sketchpad: sketching as a visual chain of thought for multimodal language models\.External Links:2406\.09403,[Link](https://arxiv.org/abs/2406.09403)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px2.p1.1)\.
- \[21\]X\. Huang, W\. Liu, X\. Zeng, Y\. Huang, X\. Hao, Y\. Wang, Y\. Zeng, C\. Wu, Y\. Wang, R\. Tang, and D\. Lian\(2025\)ToolACE\-dev: self\-improving tool learning via decomposition and evolution\.External Links:2505\.07512,[Link](https://arxiv.org/abs/2505.07512)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px3.p1.1)\.
- \[22\]W\. Ji, W\. Yuan, E\. Getzen, K\. Cho, M\. I\. Jordan, S\. Mei, J\. E\. Weston, W\. J\. Su, J\. Xu, and L\. Zhang\(2025\)An overview of large language models for statisticians\.External Links:2502\.17814,[Link](https://arxiv.org/abs/2502.17814)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px1.p1.1)\.
- \[23\]S\. A\. Joseph, S\. M\. Husain, S\. S\. Offner, S\. Juneau, P\. Torrey, A\. S\. Bolton, J\. P\. Farias, N\. Gaffney, G\. Durrett, and J\. J\. Li\(2025\)Astrovisbench: a code benchmark for scientific computing and visualization in astronomy\.arXiv preprint arXiv:2505\.20538\.Cited by:[§1](https://arxiv.org/html/2606.00384#S1.p1.1)\.
- \[24\]L\. Jung\-Mok, N\. Hyeon\-Woo, M\. Ye\-Bin, J\. Nam, and T\. Oh\(2025\)Automated model discovery via multi\-modal & multi\-step pipeline\.External Links:2509\.25946,[Link](https://arxiv.org/abs/2509.25946)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px2.p1.1)\.
- \[25\]P\. Kroupa\(2001\-04\)On the variation of the initial mass function\.Monthly Notices of the Royal Astronomical Society322\(2\),pp\. 231–246\.External Links:ISSN 1365\-2966,[Link](http://dx.doi.org/10.1046/j.1365-8711.2001.04022.x),[Document](https://dx.doi.org/10.1046/j.1365-8711.2001.04022.x)Cited by:[§4\.1](https://arxiv.org/html/2606.00384#S4.SS1.p3.1)\.
- \[26\]H\. Li, Z\. Wang, Q\. Dai, Y\. Nie, J\. Peng, R\. Liu, J\. Zhang, K\. Zhu, J\. He, L\. Wang, Y\. Ding, Y\. Chen, W\. Guo, and D\. Song\(2026\)OpenSage: self\-programming agent generation engine\.External Links:2602\.16891,[Link](https://arxiv.org/abs/2602.16891)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px3.p1.1)\.
- \[27\]M\. Y\. Li, E\. B\. Fox, and N\. D\. Goodman\(2024\)Automated statistical model discovery with language models\.External Links:2402\.17879,[Link](https://arxiv.org/abs/2402.17879)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px1.p1.1),[§4\.2](https://arxiv.org/html/2606.00384#S4.SS2.SSS0.Px1.p1.5),[§5\.2](https://arxiv.org/html/2606.00384#S5.SS2.SSS0.Px1.p1.3)\.
- \[28\]M\. Y\. Li, V\. Vajipey, N\. D\. Goodman, and E\. B\. Fox\(2024\)CriticAL: critic automation with language models\.External Links:2411\.06590,[Link](https://arxiv.org/abs/2411.06590)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px1.p1.1)\.
- \[29\]J\. R\. Lloyd, D\. Duvenaud, R\. Grosse, J\. B\. Tenenbaum, and Z\. Ghahramani\(2014\)Automatic construction and natural\-language description of nonparametric regression models\.External Links:1402\.4304,[Link](https://arxiv.org/abs/1402.4304)Cited by:[§1](https://arxiv.org/html/2606.00384#S1.p4.3)\.
- \[30\]J\. Lu, Z\. Kong, Y\. Wang, R\. Fu, H\. Wan, C\. Yang, W\. Lou, H\. Sun, L\. Wang, Y\. Jiang, X\. Wang, X\. Sun, and D\. Zhou\(2026\)Beyond static tools: test\-time tool evolution for scientific reasoning\.External Links:2601\.07641,[Link](https://arxiv.org/abs/2601.07641)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px3.p1.1)\.
- \[31\]R\. Maller, S\. Resnick, S\. Shemehsavar, and M\. Zhao\(2024\-01\)Mixture cure model methodology in survival analysis: some recent results for the one\-sample case\.Statistics Surveys18,pp\.\.External Links:[Document](https://dx.doi.org/10.1214/24-SS147)Cited by:[§4\.1](https://arxiv.org/html/2606.00384#S4.SS1.p1.3)\.
- \[32\]NASA ScienceDawn mission overview\.Note:[https://science\.nasa\.gov/mission/dawn/](https://science.nasa.gov/mission/dawn/)Accessed: May 2, 2026Cited by:[footnote 2](https://arxiv.org/html/2606.00384#footnote2)\.
- \[33\]NASA ScienceVesta: in depth\.Note:[https://science\.nasa\.gov/solar\-system/asteroids/4\-vesta/](https://science.nasa.gov/solar-system/asteroids/4-vesta/)Accessed: May 2, 2026Cited by:[footnote 1](https://arxiv.org/html/2606.00384#footnote1)\.
- \[34\]J\. Ni, Z\. Zhao, C\. Shen, H\. Tong, D\. Song, W\. Cheng, D\. Luo, and H\. Chen\(2025\)Harnessing vision models for time series analysis: a survey\.External Links:2502\.08869,[Link](https://arxiv.org/abs/2502.08869)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px1.p1.1)\.
- \[35\]S\. S\. R\. Offner, P\. C\. Clark, P\. Hennebelle, N\. Bastian, M\. R\. Bate, P\. F\. Hopkins, E\. Moreaux, and A\. P\. Whitworth\(2014\)The origin and universality of the stellar initial mass function\.InProtostars and Planets VI,External Links:ISBN 9780816531240,[Link](http://dx.doi.org/10.2458/azu_uapress_9780816531240-ch003),[Document](https://dx.doi.org/10.2458/azu%5Fuapress%5F9780816531240-ch003)Cited by:[§4\.1](https://arxiv.org/html/2606.00384#S4.SS1.p2.1)\.
- \[36\]C\. Qian, C\. Han, Y\. R\. Fung, Y\. Qin, Z\. Liu, and H\. Ji\(2024\)CREATOR: tool creation for disentangling abstract and concrete reasoning of large language models\.External Links:2305\.14318,[Link](https://arxiv.org/abs/2305.14318)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px3.p1.1)\.
- \[37\]J\. Requeima, J\. Bronskill, D\. Choi, R\. E\. Turner, and D\. Duvenaud\(2024\)LLM processes: numerical predictive distributions conditioned on natural language\.External Links:2405\.12856,[Link](https://arxiv.org/abs/2405.12856)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px1.p1.1)\.
- \[38\]W\. Rudman, M\. Golovanevsky, A\. Bar, V\. Palit, Y\. LeCun, C\. Eickhoff, and R\. Singh\(2025\-07\)Forgotten polygons: multimodal large language models are shape\-blind\.InFindings of the Association for Computational Linguistics: ACL 2025,W\. Che, J\. Nabende, E\. Shutova, and M\. T\. Pilehvar \(Eds\.\),Vienna, Austria,pp\. 11983–11998\.External Links:[Link](https://aclanthology.org/2025.findings-acl.620/),[Document](https://dx.doi.org/10.18653/v1/2025.findings-acl.620),ISBN 979\-8\-89176\-256\-5Cited by:[§6\.2](https://arxiv.org/html/2606.00384#S6.SS2.SSS0.Px1.p1.1)\.
- \[39\]E\. E\. Salpeter\(1955\-01\)The Luminosity Function and Stellar Evolution\.\.apj121,pp\. 161\.External Links:[Document](https://dx.doi.org/10.1086/145971)Cited by:[§4\.1](https://arxiv.org/html/2606.00384#S4.SS1.p3.1)\.
- \[40\]C\. Si, Z\. Yang, Y\. Choi, E\. Candès, D\. Yang, and T\. Hashimoto\(2026\)Towards execution\-grounded automated ai research\.External Links:2601\.14525,[Link](https://arxiv.org/abs/2601.14525)Cited by:[§1](https://arxiv.org/html/2606.00384#S1.p2.1)\.
- \[41\]Y\. Song, W\. Xiong, D\. Zhu, W\. Wu, H\. Qian, M\. Song, H\. Huang, C\. Li, K\. Wang, R\. Yao, Y\. Tian, and S\. Li\(2023\)RestGPT: connecting large language models with real\-world restful apis\.External Links:2306\.06624,[Link](https://arxiv.org/abs/2306.06624)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px3.p1.1)\.
- \[42\]M\. Sun, R\. Han, B\. Jiang, H\. Qi, D\. Sun, Y\. Yuan, and J\. Huang\(2025\-10\)A survey on large language model\-based agents for statistics and data science\.The American Statistician,pp\. 1–14\.External Links:ISSN 1537\-2731,[Link](http://dx.doi.org/10.1080/00031305.2025.2561140),[Document](https://dx.doi.org/10.1080/00031305.2025.2561140)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px1.p1.1)\.
- \[43\]Z\. Sun, Z\. Liu, Y\. Zang, Y\. Cao, X\. Dong, T\. Wu, D\. Lin, and J\. Wang\(2025\)SEAgent: self\-evolving computer use agent with autonomous learning from experience\.External Links:2508\.04700,[Link](https://arxiv.org/abs/2508.04700)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px3.p1.1)\.
- \[44\]D\. Surís, S\. Menon, and C\. Vondrick\(2023\)ViperGPT: visual inference via python execution for reasoning\.External Links:2303\.08128,[Link](https://arxiv.org/abs/2303.08128)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px2.p1.1)\.
- \[45\]S\. Tong, Z\. Liu, Y\. Zhai, Y\. Ma, Y\. LeCun, and S\. Xie\(2024\)Eyes wide shut? exploring the visual shortcomings of multimodal llms\.External Links:2401\.06209,[Link](https://arxiv.org/abs/2401.06209)Cited by:[§6\.2](https://arxiv.org/html/2606.00384#S6.SS2.SSS0.Px1.p1.1)\.
- \[46\]G\. Wang, Y\. Xie, Y\. Jiang, A\. Mandlekar, C\. Xiao, Y\. Zhu, L\. Fan, and A\. Anandkumar\(2023\)Voyager: an open\-ended embodied agent with large language models\.External Links:2305\.16291,[Link](https://arxiv.org/abs/2305.16291)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px3.p1.1)\.
- \[47\]Q\. Wen, T\. Zhou, C\. Zhang, W\. Chen, Z\. Ma, J\. Yan, and L\. Sun\(2023\)Transformers in time series: a survey\.External Links:2202\.07125,[Link](https://arxiv.org/abs/2202.07125)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px1.p1.1)\.
- \[48\]G\. Wölflein, D\. Ferber, D\. Truhn, O\. Arandjelović, and J\. N\. Kather\(2025\)LLM agents making agent tools\.External Links:2502\.11705,[Link](https://arxiv.org/abs/2502.11705)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px3.p1.1)\.
- \[49\]C\. Wu, S\. Yin, W\. Qi, X\. Wang, Z\. Tang, and N\. Duan\(2023\)Visual chatgpt: talking, drawing and editing with visual foundation models\.External Links:2303\.04671,[Link](https://arxiv.org/abs/2303.04671)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px2.p1.1)\.
- \[50\]S\. Yan, J\. Tong, H\. Xue, X\. Tang, Y\. Wang, K\. Shi, G\. Zhang, R\. Li, and Y\. Zou\(2026\)Act wisely: cultivating meta\-cognitive tool use in agentic multimodal models\.External Links:2604\.08545,[Link](https://arxiv.org/abs/2604.08545)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px2.p1.1)\.
- \[51\]X\. Yu, C\. Xu, G\. Zhang, Z\. Chen, Y\. Zhang, Y\. He, P\. Jiang, J\. Zhang, X\. Hu, and S\. Yan\(2026\)VisMem: latent vision memory unlocks potential of vision\-language models\.External Links:2511\.11007,[Link](https://arxiv.org/abs/2511.11007)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px2.p1.1)\.
- \[52\]G\. Zerveas, S\. Jayaraman, D\. Patel, A\. Bhamidipaty, and C\. Eickhoff\(2020\)A transformer\-based framework for multivariate time series representation learning\.External Links:2010\.02803,[Link](https://arxiv.org/abs/2010.02803)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px1.p1.1)\.
- \[53\]Y\. Zhang, L\. Hu, H\. Sun, P\. Wang, Y\. Wei, S\. Yin, J\. Pei, W\. Shen, P\. Xia, Y\. Peng, T\. Xie, E\. Li, Y\. Liu, X\. Song, and Y\. Zhou\(2025\)Skywork\-r1v4: toward agentic multimodal intelligence through interleaved thinking with images and deepresearch\.External Links:2512\.02395,[Link](https://arxiv.org/abs/2512.02395)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px2.p1.1)\.
- \[54\]Z\. Zhang, R\. Rossi, T\. Yu, F\. Dernoncourt, R\. Zhang, J\. Gu, S\. Kim, X\. Chen, Z\. Wang, and N\. Lipka\(2025\)VipAct: visual\-perception enhancement via specialized vlm agent collaboration and tool\-use\.External Links:2410\.16400,[Link](https://arxiv.org/abs/2410.16400)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px2.p1.1)\.
- \[55\]S\. Zhao, H\. Zhang, S\. Lin, M\. Li, Q\. Wu, K\. Zhang, and C\. Wei\(2025\)PyVision: agentic vision with dynamic tooling\.External Links:2507\.07998,[Link](https://arxiv.org/abs/2507.07998)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px2.p1.1),[§5\.2](https://arxiv.org/html/2606.00384#S5.SS2.SSS0.Px2.p1.1)\.
- \[56\]S\. Zhong, W\. Ruan, M\. Jin, H\. Li, Q\. Wen, and Y\. Liang\(2025\)Time\-vlm: exploring multimodal vision\-language models for augmented time series forecasting\.External Links:2502\.04395,[Link](https://arxiv.org/abs/2502.04395)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px1.p1.1)\.
- \[57\]Q\. Zhou, R\. Zhou, Z\. Hu, P\. Lu, S\. Gao, and Y\. Zhang\(2024\)Image\-of\-thought prompting for visual reasoning refinement in multimodal large language models\.External Links:2405\.13872,[Link](https://arxiv.org/abs/2405.13872)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px2.p1.1)\.
- \[58\]Z\. Zhou, D\. Chen, Z\. Ma, Z\. Hu, M\. Fu, S\. Wang, Y\. Wan, Z\. Zhao, and R\. Krishna\(2025\)Reinforced visual perception with tools\.External Links:2509\.01656,[Link](https://arxiv.org/abs/2509.01656)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px2.p1.1)\.
- \[59\]Y\. Zhu, S\. Du, B\. Li, Y\. Luo, and N\. Tang\(2024\)Are large language models good statisticians?\.External Links:2406\.07815,[Link](https://arxiv.org/abs/2406.07815)Cited by:[§2](https://arxiv.org/html/2606.00384#S2.SS0.SSS0.Px1.p1.1)\.

## Appendix AVestaDetails

Algorithm 2Visual Exploration Agents \(Detailed\)1:Data

DD, iteration limit

NN, proposals per iteration

pp, metric

RR, registry

ℰ\\mathcal\{E\}\(initial state:generate\_new\_toolonly\)

2:

Mbest,θbestM\_\{\\text\{best\}\},\\theta\_\{\\text\{best\}\}
3:

P0←Plot​\(D\)P\_\{0\}\\leftarrow\\text\{Plot\}\(D\)
4:

\{M0\(j\)\}j=1p←Propose​\(P0\)\\\{M\_\{0\}^\{\(j\)\}\\\}\_\{j=1\}^\{p\}\\leftarrow\\textsc\{Propose\}\(P\_\{0\}\)⊳\\trianglerightProposal Agent proposesppdiverse models

5:

θ0\(j\)←StatGen​\(D,M0\(j\)\)​∀j\\theta\_\{0\}^\{\(j\)\}\\leftarrow\\textsc\{StatGen\}\(D,M\_\{0\}^\{\(j\)\}\)\\ \\forall j⊳\\trianglerightStat\-Gen fits allppmodels

6:

M0,θ0←arg⁡minj⁡R​\(D,M0\(j\)​\(θ0\(j\)\)\)M\_\{0\},\\theta\_\{0\}\\leftarrow\\arg\\min\_\{j\}\\ R\(D,M\_\{0\}^\{\(j\)\}\(\\theta\_\{0\}^\{\(j\)\}\)\)⊳\\trianglerightSelect best model under metricRR

7:

Mbest,θbest←M0,θ0M\_\{\\text\{best\}\},\\theta\_\{\\text\{best\}\}\\leftarrow M\_\{0\},\\theta\_\{0\}
8:

s0←Summarize​\(Mbest,θbest\)s\_\{0\}\\leftarrow\\textsc\{Summarize\}\(M\_\{\\text\{best\}\},\\theta\_\{\\text\{best\}\}\)
9:for

i=1,…,Ni=1,\\ldots,Ndo

10:

ε∗←ToolManager​\(si−1,D,Mi−1,θi−1,ℰ\)\\varepsilon^\{\*\}\\leftarrow\\textsc\{ToolManager\}\(s\_\{i\-1\},D,M\_\{i\-1\},\\theta\_\{i\-1\},\\mathcal\{E\}\)⊳\\trianglerightSelect orcreatediagnostic tool

11:

ℰ←ℰ∪\{ε∗\}\\mathcal\{E\}\\leftarrow\\mathcal\{E\}\\cup\\\{\\varepsilon^\{\*\}\\\}⊳\\trianglerightPersist newly created tools to registry

12:

Ti←ε∗​\(D,Mi−1,θi−1\)T\_\{i\}\\leftarrow\\varepsilon^\{\*\}\(D,M\_\{i\-1\},\\theta\_\{i\-1\}\)⊳\\trianglerightExecute tool to produce diagnostic

13:

desc,\{Mi\(j\)\}j=1p←Critique​\(Ti,s0:i−1,D\)\\text\{desc\},\\ \\\{M\_\{i\}^\{\(j\)\}\\\}\_\{j=1\}^\{p\}\\leftarrow\\textsc\{Critique\}\(T\_\{i\},s\_\{0:i\-1\},D\)⊳\\trianglerightProposepprevised models

14:if

desc=COMPLETE\\text\{desc\}=\\texttt\{COMPLETE\}then

15:break

16:endif

17:

θi\(j\)←StatGen​\(D,Mi\(j\)\)​∀j\\theta\_\{i\}^\{\(j\)\}\\leftarrow\\textsc\{StatGen\}\(D,M\_\{i\}^\{\(j\)\}\)\\ \\forall j⊳\\trianglerightStat\-Gen fits allppmodels

18:

Mi,θi←arg⁡minj⁡R​\(D,Mi\(j\)​\(θi\(j\)\)\)M\_\{i\},\\theta\_\{i\}\\leftarrow\\arg\\min\_\{j\}\\ R\(D,M\_\{i\}^\{\(j\)\}\(\\theta\_\{i\}^\{\(j\)\}\)\)⊳\\trianglerightSelect best model under metricRR

19:

Mbest,θbest←arg⁡minθ∈\{θi,θbest\}⁡R​\(D,ℳ​\(θ\)\)M\_\{\\text\{best\}\},\\theta\_\{\\text\{best\}\}\\leftarrow\\arg\\min\_\{\\theta\\in\\\{\\theta\_\{i\},\\,\\theta\_\{\\text\{best\}\}\\\}\}R\(D,\\mathcal\{M\}\(\\theta\)\)
20:

si←Summarize​\(Mi,θi,ε∗\)s\_\{i\}\\leftarrow\\textsc\{Summarize\}\(M\_\{i\},\\theta\_\{i\},\\varepsilon^\{\*\}\)
21:endfor

22:return

Mbest,θbestM\_\{\\text\{best\}\},\\theta\_\{\\text\{best\}\}

## Appendix BDawnDetails\.

### B\.1Distribution Fitting

For eachDTrueiD\_\{\\text\{True\}\_\{i\}\}in the distribution fitting domain, ground\-truth parameters are drawn independently and uniformly from predefined ranges, with the family \(or mixture of families\) depending on the difficulty split\. Eight base families are used throughout: Gaussian, Log\-normal, Student\-tt, Exponential, Uniform, Weibull, Cauchy, and Laplace\.

##### Easy split\.

A single distribution family is drawn uniformly at random from the eight families above\. Parameters are then sampled as follows\. For the*Gaussian*, meanμ∼𝒰​\(−5,5\)\\mu\\sim\\mathcal\{U\}\(\-5,5\)and standard deviationσ∼𝒰​\(0\.5,3\)\\sigma\\sim\\mathcal\{U\}\(0\.5,3\)\. For the*Log\-normal*, log\-meanμ∼𝒰​\(0,2\)\\mu\\sim\\mathcal\{U\}\(0,2\)and log\-widthσ∼𝒰​\(0\.3,1\.5\)\\sigma\\sim\\mathcal\{U\}\(0\.3,1\.5\)\. For the*Student\-tt*, locationμ∼𝒰​\(−5,5\)\\mu\\sim\\mathcal\{U\}\(\-5,5\), scaleσ∼𝒰​\(0\.5,3\)\\sigma\\sim\\mathcal\{U\}\(0\.5,3\), and degrees of freedomν∼𝒰​\(2\.5,10\)\\nu\\sim\\mathcal\{U\}\(2\.5,10\)\. For the*Exponential*, locationℓ∼𝒰​\(0,2\)\\ell\\sim\\mathcal\{U\}\(0,2\)and scaleλ∼𝒰​\(0\.5,3\)\\lambda\\sim\\mathcal\{U\}\(0\.5,3\)\. For the*Uniform*, meanμ∼𝒰​\(−5,5\)\\mu\\sim\\mathcal\{U\}\(\-5,5\)and standard deviationσ∼𝒰​\(0\.5,3\)\\sigma\\sim\\mathcal\{U\}\(0\.5,3\), from which the bounds are derived as\[μ−σ​3,μ\+σ​3\]\[\\mu\-\\sigma\\sqrt\{3\},\\;\\mu\+\\sigma\\sqrt\{3\}\]\. For the*Weibull*, locationℓ∼𝒰​\(0,2\)\\ell\\sim\\mathcal\{U\}\(0,2\), shapeα∼𝒰​\(1,3\)\\alpha\\sim\\mathcal\{U\}\(1,3\), and scaleλ∼𝒰​\(0\.5,2\)\\lambda\\sim\\mathcal\{U\}\(0\.5,2\)\. For the*Cauchy*, locationα∼𝒰​\(−5,5\)\\alpha\\sim\\mathcal\{U\}\(\-5,5\)and scaleβ∼𝒰​\(0\.5,2\)\\beta\\sim\\mathcal\{U\}\(0\.5,2\)\. For the*Laplace*, locationμ∼𝒰​\(−5,5\)\\mu\\sim\\mathcal\{U\}\(\-5,5\)and diversityb∼𝒰​\(0\.5,2\)b\\sim\\mathcal\{U\}\(0\.5,2\)\. Each dataset containsn∈\[500,1500\]n\\in\[500,1500\]points\.

##### Hard split\.

Two distribution families are drawn independently with replacement from the same eight families\. Each component draws its own sample sizenk∈\[250,750\]n\_\{k\}\\in\[250,750\]independently, and parameters for each component are sampled using the same per\-family ranges described above\. The two components are concatenated to form the final dataset of total sizen=n1\+n2n=n\_\{1\}\+n\_\{2\}\.

##### Astro split\.

Each ground\-truth model is one of five initial mass function \(IMF\) families: Salpeter, Kroupa, Chabrier, Free\-form \(wide\), and Free\-form \(tight\), with 20 independent realisations generated per family\. The number of stellar masses drawn per realisation is sampled uniformly from\[500,1500\]\[500,1500\], and all families sharemmax=120​M⊙m\_\{\\max\}=120\\,M\_\{\\odot\}\. The Salpeter IMF is a single power lawξ​\(m\)∝m−α\\xi\(m\)\\propto m^\{\-\\alpha\}withα∼𝒰​\(1\.9,2\.7\)\\alpha\\sim\\mathcal\{U\}\(1\.9,2\.7\)on\[0\.08,120\]​M⊙\[0\.08,120\]\\,M\_\{\\odot\}\. The Kroupa IMF is a three\-segment broken power law with slopesp1∼𝒰​\(0\.1,0\.5\)p\_\{1\}\\sim\\mathcal\{U\}\(0\.1,0\.5\),p2∼𝒰​\(1\.0,1\.6\)p\_\{2\}\\sim\\mathcal\{U\}\(1\.0,1\.6\),p3∼𝒰​\(1\.9,2\.7\)p\_\{3\}\\sim\\mathcal\{U\}\(1\.9,2\.7\)and fixed break massesm1=0\.08​M⊙m\_\{1\}=0\.08\\,M\_\{\\odot\},m2=0\.5​M⊙m\_\{2\}=0\.5\\,M\_\{\\odot\}on\[0\.009,120\]​M⊙\[0\.009,120\]\\,M\_\{\\odot\}; continuity constants are enforced at each break\. The Chabrier IMF joins a log\-normal component \(mchar∼𝒰​\(0\.15,0\.35\)m\_\{\\text\{char\}\}\\sim\\mathcal\{U\}\(0\.15,0\.35\),σ∼𝒰​\(0\.45,0\.70\)\\sigma\\sim\\mathcal\{U\}\(0\.45,0\.70\)\) below1​M⊙1\\,M\_\{\\odot\}to a power law \(α∼𝒰​\(1\.9,2\.7\)\\alpha\\sim\\mathcal\{U\}\(1\.9,2\.7\)\) above1​M⊙1\\,M\_\{\\odot\}, matched continuously at the junction\. The Free\-form IMF is an 8\-parameter broken power law with four slopes and three interior break masses; the wide variant uses broad exploratory ranges \(a0∼𝒰​\(−0\.5,1\.0\)a\_\{0\}\\\!\\sim\\\!\\mathcal\{U\}\(\-0\.5,1\.0\),a1∼𝒰​\(0\.5,2\.0\)a\_\{1\}\\\!\\sim\\\!\\mathcal\{U\}\(0\.5,2\.0\),a2∼𝒰​\(1\.5,3\.0\)a\_\{2\}\\\!\\sim\\\!\\mathcal\{U\}\(1\.5,3\.0\),a3∼𝒰​\(1\.8,3\.5\)a\_\{3\}\\\!\\sim\\\!\\mathcal\{U\}\(1\.8,3\.5\); breaksb1∼𝒰​\(0\.06,0\.12\)b\_\{1\}\\\!\\sim\\\!\\mathcal\{U\}\(0\.06,0\.12\),b2∼𝒰​\(0\.30,0\.80\)b\_\{2\}\\\!\\sim\\\!\\mathcal\{U\}\(0\.30,0\.80\),b3∼𝒰​\(0\.80,2\.00\)b\_\{3\}\\\!\\sim\\\!\\mathcal\{U\}\(0\.80,2\.00\);mmin∼𝒰​\(0\.007,0\.040\)m\_\{\\min\}\\\!\\sim\\\!\\mathcal\{U\}\(0\.007,0\.040\)\), while the tight variant uses ranges anchored near canonical Kroupa values \(a0∼𝒰​\(0\.1,0\.5\)a\_\{0\}\\\!\\sim\\\!\\mathcal\{U\}\(0\.1,0\.5\),a1∼𝒰​\(1\.0,1\.6\)a\_\{1\}\\\!\\sim\\\!\\mathcal\{U\}\(1\.0,1\.6\),a2,a3∼𝒰​\(1\.9,2\.7\)a\_\{2\},a\_\{3\}\\\!\\sim\\\!\\mathcal\{U\}\(1\.9,2\.7\); breaksb1∼𝒰​\(0\.07,0\.09\)b\_\{1\}\\\!\\sim\\\!\\mathcal\{U\}\(0\.07,0\.09\),b2∼𝒰​\(0\.40,0\.60\)b\_\{2\}\\\!\\sim\\\!\\mathcal\{U\}\(0\.40,0\.60\),b3∼𝒰​\(0\.80,1\.20\)b\_\{3\}\\\!\\sim\\\!\\mathcal\{U\}\(0\.80,1\.20\);mmin∼𝒰​\(0\.007,0\.012\)m\_\{\\min\}\\\!\\sim\\\!\\mathcal\{U\}\(0\.007,0\.012\)\)\. Break masses are sorted to enforceb1<b2<b3b\_\{1\}<b\_\{2\}<b\_\{3\}, and continuity constants are applied at each break\.

### B\.2Time Series

For eachDTrueiD\_\{\\text\{True\}\_\{i\}\}in the time\-series domain, ground\-truth signal parameters are drawn independently and uniformly from predefined ranges, with the specific ranges and signal family depending on the difficulty split\.

##### Easy split\.

Each ground\-truth signal follows a sinusoidal trend model,

y​\(t\)=s⋅t~\+A​sin⁡\(2​π​f​t~\)\+ε​\(t\),y\(t\)=s\\cdot\\tilde\{t\}\+A\\sin\(2\\pi f\\tilde\{t\}\)\+\\varepsilon\(t\),\(2\)wheret~∈\[0,1\]\\tilde\{t\}\\in\[0,1\]is normalized time\. Parameters are sampled as slopes∼𝒰​\(−2\.5,2\.5\)s\\sim\\mathcal\{U\}\(\-2\.5,\\,2\.5\), amplitudeA∼𝒰​\(0\.2,2\.0\)A\\sim\\mathcal\{U\}\(0\.2,\\,2\.0\), frequencyf∼𝒰​\(1,12\)f\\sim\\mathcal\{U\}\(1,\\,12\)cycles per window, and noise fractionσ/A∼𝒰​\(0\.05,0\.15\)\\nicefrac\{\{\\sigma\}\}\{\{A\}\}\\sim\\mathcal\{U\}\(0\.05,\\,0\.15\)withε​\(t\)∼𝒩​\(0,σ2\)\\varepsilon\(t\)\\sim\\mathcal\{N\}\(0,\\sigma^\{2\}\)\. Approximately 16% of series are further constrained tos=0s=0to ensure flat\-trend instances are represented\. All series are mean\-centred after generation\.

##### Hard split\.

A two\-level sampling scheme is used\. First, a signal category is drawn—Complexwith probability0\.40\.4orRhythmicwith probability0\.60\.6\. WithinComplex, a type is drawn uniformly from: ARMA \(an AR\(1\) random walk\), Heteroscedastic \(an expanding\-variance process with additive seasonal component\), and Logistic \(an S\-curve trend\)\. WithinRhythmic, a type is drawn uniformly from: synthetic ECG \(a Gaussian P\-QRS\-T complex template\) and Square wave, with period sampled from\{40,60,80\}\\\{40,60,80\\\}time steps; small Gaussian noiseε​\(t\)∼𝒩​\(0,σ2\)\\varepsilon\(t\)\\sim\\mathcal\{N\}\(0,\\sigma^\{2\}\)is added to all rhythmic signals\. Pure sine\-wave instances are generated but excluded from the final dataset to prevent overlap with the easy split\.

##### Astro split \(gravitational wave chirps\)\.

Each ground\-truth signal follows a linearly swept\-frequency \(chirp\) model,

y​\(t\)=A​\(t\)​sin⁡\(2​π​\(f0​t\+12​β​t2\)\)\+ε​\(t\),y\(t\)=A\(t\)\\sin\\\!\\left\(2\\pi\\\!\\left\(f\_\{0\}\\,t\+\\tfrac\{1\}\{2\}\\beta\\,t^\{2\}\\right\)\\right\)\+\\varepsilon\(t\),\(3\)giving instantaneous frequencyf​\(t\)=f0\+β​tf\(t\)=f\_\{0\}\+\\beta t\. Parameters are sampled asf0∼𝒰​\(0\.3,2\.5\)f\_\{0\}\\sim\\mathcal\{U\}\(0\.3,\\,2\.5\),β∼𝒰​\(0\.20,1\.00\)\\beta\\sim\\mathcal\{U\}\(0\.20,\\,1\.00\),A∼𝒰​\(0\.7,1\.8\)A\\sim\\mathcal\{U\}\(0\.7,\\,1\.8\), andσ∼𝒰​\(0\.07,0\.15\)\\sigma\\sim\\mathcal\{U\}\(0\.07,\\,0\.15\)\. For 40 of the 50 series, amplitude is constant,A​\(t\)=AA\(t\)=A; for the remaining 10, an exponential decay envelopeA​\(t\)=A​e−t/τA\(t\)=A\\,e^\{\-t/\\tau\}withτ∼𝒰​\(2\.0,5\.0\)\\tau\\sim\\mathcal\{U\}\(2\.0,\\,5\.0\)is applied to model inspiral ringdown behaviour\.

All time series consist ofN=600N=600equally spaced observations on a daily grid\.

## Appendix CVestaTool Generation Variants

To evaluate whether tools can be effectively accumulated through sequential interaction, we sampled 100 instances from theHarddistribution fitting split, partitioning them equally into 50 train and 50 test examples\. We evaluate three conditions with Claude Sonnet 4\.6\. \(1\) Accumulate\-Only: During a train phase, the model processes training examples sequentially with tool accumulation enabled, initializing each run with the tools generated in all prior runs and selecting either an existing tool or generating a new one\. After 50 sequential runs, this process converges to just 4 unique tools\. During the test phase, the model is restricted to selecting from these accumulated tools with no new generation permitted and no interaction between runs\. \(2\) Generate\-Only: Identical to Accumulate\-Only, but the test phase additionally permits the model to generate new tools, which are not shared across test runs\. \(3\) Expert Toolkit \+ Generation: No train phase is used\. Instead, the model is initialized with the default expert toolkit on all 50 test examples while also retaining the ability to generate new tools, again without cross\-run sharing\. This design isolates three factors: the value of accumulated tools alone, the benefit of combining accumulated tools with on\-the\-fly generation, and the contribution of generation when paired with a fixed expert toolkit\.

![Refer to caption](https://arxiv.org/html/2606.00384v1/images/accumulated_variant_comparison.png)Figure 6:Jensen–Shannon divergence on theHarddistribution fitting split \(lower is better\) comparing three accumulated toolkit conditions with Claude Sonnet 4\.6\.presents JS divergence scores across the three toolkit conditions on the Hard distribution fitting split\. Overall, performance differences across conditions are modest, with all three variants achieving mean JS divergence between 0\.106 and 0\.124, and overlapping standard errors indicating no statistically significant separation\. Nevertheless, a consistent ordering emerges: Expert\+Generate achieves the lowest divergence \(0\.106\), followed by Generate\-Only \(0\.115\), with Accumulate\-Only performing worst \(0\.124\)\. Two conclusions follow from this pattern\. First, tool accumulation alone provides no measurable benefit\. Restricting the model to the four auto\-accumulated tools at test time yields the weakest performance, suggesting that sequential accumulation over 50 training examples converges to a toolkit that is too narrow to generalize across the full diversity of dataset\. Second, on\-the\-fly tool generation is the primary performance driver\. Both conditions that permit test\-time generation outperform Accumulate\-Only, regardless of whether the base toolkit is accumulated or expert\-designed\. The marginal advantage of Expert\+Generate over Generate\-Only further suggests that a well\-curated expert toolkit provides a stronger generative prior than one derived through automatic accumulation\. Taken together, these results indicate that the current accumulation strategy does not yet yield tools competitive with expert\-designed ones, and that allowing the model to generate new tools at test time remains essential for robust distribution fitting performance\.

## Appendix DAdditional Results onVesta

### D\.1Statistical Significance of the Results

We pair\-wise test, for each combination of LLM, dataset split, baseline, andVestavariant, whetherVestasignificantly outperforms the baseline\. The two evaluation metrics are treated as independent statistical questions: they probe different tasks \(distribution fitting vs\. time\-series modelling\), have different data structures, and warrant different tests\. Multiple\-comparison corrections are applied within each metric but never across them\.

##### Distribution fitting \(JS divergence\)\.

Per\-problem JS values are available for every method on the same set of problems, so we use a one\-sided paired Wilcoxon signed\-rank test across problems withH1:JSVesta<JSbaselineH\_\{1\}\\\!:\\mathrm\{JS\}\_\{\{\\sc Vesta\}\}<\\mathrm\{JS\}\_\{\\text\{baseline\}\}\.

##### Time\-series modeling \(ELPD LOO\)\.

ArviZ’slooreturns a point estimateELPD^\\widehat\{\\mathrm\{ELPD\}\}and its standard errorSE^\\widehat\{\\mathrm\{SE\}\}for each fit\. We use a one\-sidedzz\-test onΔ​ELPD=ELPDVesta−ELPDbaseline\\Delta\\mathrm\{ELPD\}\\\!=\\\!\\mathrm\{ELPD\}\_\{\{\\sc Vesta\}\}\\\!\-\\\!\\mathrm\{ELPD\}\_\{\\text\{baseline\}\}withSEΔ=SE^V2\+SE^B2\\mathrm\{SE\}\_\{\\Delta\}=\\sqrt\{\\widehat\{\\mathrm\{SE\}\}\_\{V\}^\{\\,2\}\+\\widehat\{\\mathrm\{SE\}\}\_\{B\}^\{\\,2\}\}\. ThisSEΔ\\mathrm\{SE\}\_\{\\Delta\}is*conservative*:loo\_comparecomputes the SE from per\-observation ELPD differences, which is strictly smaller whenever the two ELPDs are positively correlated\. Our reportedpp\-values are therefore upper bounds on the truepp\-values\.

##### Multiple comparisons\.

Within each metric we conduct3​splits×3​LLMs×3​Vestavariants×2​baselines=543\\text\{ splits\}\\times 3\\text\{ LLMs\}\\times 3\\text\{ \{\\sc Vesta\}\{\} variants\}\\times 2\\text\{ baselines\}=54tests and apply Holm–Bonferroni correction\. We useα=0\.05\\alpha=0\.05throughout\.

### D\.2Summary: aggregate win rates

Table[4](https://arxiv.org/html/2606.00384#A4.T4)reports, per \(variant, baseline, split\), how many of the 3 LLMs yield a statistically significant win forVestaafter Holm correction\.Vesta’s advantage on time\-series modelling is robust: every variant beats PyVision in every \(LLM, split\) cell, and beats Box\-LM in 13/18 cells\. For distribution fitting, significant wins are concentrated on theHardsplit, where the absolute JS gap betweenVestaand the baselines is largest\.

Table 4:Number of LLMs \(out of 3\) for which eachVestavariant significantly outperforms each baseline atα=0\.05\\alpha\\\!=\\\!0\.05after Holm correction\. Tests are one\-sided in the direction “Vestabeats baseline” \(paired Wilcoxon for JS divergence;zz\-test onΔ\\DeltaELPD for ELPD LOO\)\. Bold marks cells with≥1\\geq\\\!1significant win\.
### D\.3Detailed per\-LLM results

Tables[5](https://arxiv.org/html/2606.00384#A4.T5)and[6](https://arxiv.org/html/2606.00384#A4.T6)give the median JS gap andΔ\\DeltaELPD respectively for every individual comparison, along with stars denoting Holm\-corrected significance levels\.

Table 5:Distribution fitting: median JS gapΔ=med​\(JSVesta\)−med​\(JSbaseline\)\\Delta\\\!=\\\!\\mathrm\{med\}\(\\mathrm\{JS\}\_\{\{\\sc Vesta\}\}\)\-\\mathrm\{med\}\(\\mathrm\{JS\}\_\{\\text\{baseline\}\}\)\. Negative values favourVesta\. Stars are Holm\-correctedpp\-values from the one\-sided paired Wilcoxon test:\*p<0\.05p\\\!<\\\!0\.05,\*\*p<0\.01p\\\!<\\\!0\.01,\*\*\*p<0\.001p\\\!<\\\!0\.001\. Significant wins \(pHolm<0\.05p\_\{\\text\{Holm\}\}\\\!<\\\!0\.05\) in bold\.Table 6:Time\-series modelling:Δ​ELPD=ELPDVesta−ELPDbaseline\\Delta\\mathrm\{ELPD\}\\\!=\\\!\\mathrm\{ELPD\}\_\{\{\\sc Vesta\}\}\\\!\-\\\!\\mathrm\{ELPD\}\_\{\\text\{baseline\}\}, in nats \(positive favoursVesta\)\. Stars are Holm\-correctedpp\-values from the one\-sidedzz\-test:\*p<0\.05p\\\!<\\\!0\.05,\*\*p<0\.01p\\\!<\\\!0\.01,\*\*\*p<0\.001p\\\!<\\\!0\.001\. Significant wins in bold\.
### D\.4Discussion

##### Time\-series modelling\.

Vesta’s advantage on ELPD LOO is essentially universal: 45 of 54 comparisons are significant atα=0\.05\\alpha\\\!=\\\!0\.05after Holm correction, with the great majority atp<0\.001p\\\!<\\\!0\.001\. The five non\-significant comparisons against Box\-LM on Claude \(Easy and Hard splits\) are an artefact of Box\-LM’s catastrophically diffuse posterior on Claude — theΔ​ELPD\\Delta\\mathrm\{ELPD\}point estimates are massive \(\+398\+398to\+644\+644nats in favour ofVesta\), but the standard error on the Box\-LM side \(∼360\\sim\\\!360nats\) is large enough to keep these just above the Holm\-corrected threshold\. Substantively,Vestaclearly wins here; the test is simply unable to certify it given how unstable Box\-LM’s posterior is\.

##### Distribution fitting\.

Significant wins on JS divergence are concentrated in theHardsplit, where the effect size is largest\. OnEasyandAstro, point estimates favourVestain nearly every cell \(Table[5](https://arxiv.org/html/2606.00384#A4.T5)\), but per\-problem variance is high relative to the typical gap, and the Wilcoxon test cannot certify the difference atα=0\.05\\alpha\\\!=\\\!0\.05given 54 corrections\. The Kimi/Astroresults are particularly noisy:Vesta\-Dynamic posts a−0\.370\-0\.370median JS gap vs\. PyVision yetpHolm=1p\_\{\\text\{Holm\}\}\\\!=\\\!1, because the paired sample is small and skewed\. We therefore present the JS divergence results as suggestive rather than confirmatory onEasyandAstro, and as confirmatory onHard\.

Table 7:Distribution fitting: posterior quality metrics by LLM, method, and dataset\. ELPD LOO \(↑\\uparrow\): posterior predictive log\-density on held\-out data\. Bhatt \(↓\\downarrow\): Bhattacharyya distance between predicted and true distributions\.ploop\_\{\\text\{loo\}\}\(↓\\downarrow\): PSIS\-LOO effective number of parameters \(lower indicates cleaner posterior geometry\)\. Bold marks the best value within each \(LLM, dataset, metric\) column\.

### D\.5Alternative Evaluation Metrics

Table 8:Mean performance comparison across all toolkit configurations and datasets for time series\. Metrics: ELPD LOO \(higher↑\\uparrow\),R2R^\{2\}\(higher↑\\uparrow\), CRPS \(lower↓\\downarrow\); We usePyVisionas our baseline\.Estimating the Expected Log Predictive Density \(ELPD\) for Gaussian Process based Time Series models inPyMCis computationally demanding because it necessitates full MCMC sampling\. Unlike simpler parametric models, GPs require the inversion of ann×nn\\times ncovariance matrix \- a process withO​\(n3\)O\(n^\{3\}\)complexity \- at every gradient evaluation within the NUTS sampler\. To mitigate this cubic scaling bottleneck, we employ fixed\-space sub\-sampling\. By selecting a representative subset of the time series data at regular intervals, we reduce the effectivennand this significantly speeds up the evaluation\. For consistency, we measure theR2R^\{2\}, CRPS are calculated in the same manner\.

For the baseline evaluations, BoxLM generally performs better than PyVision, with the exception of a small number of anomalous cases in Claude\-generated PyMC code\. Specifically, a few extremely large ELPD values substantially inflate the mean performance metrics\. Excluding these 7 outliers \(out of 630 total samples\) reduces the mean ELPD for Claude Sonnet 4\.6 \+ BoxLM to40\.4on the Easy dataset and\-150\.6on the Hard dataset, bringing performance closer for Easy, and ahead for Hard dataset to the corresponding No\-Toolkit configuration\.

In contrast, PyVision exhibits substantially higher failure rates during model fitting\. Across datasets, approximately 40–60% of generated samples fail either due to code execution issues or unsuccessful MCMC fitting, representing the highest failure rate among all evaluated approaches\. In the remaining successful runs, the fitted models frequently collapse to near mean\-value predictions, resulting in poor overall ELPD and near\-zeroR2R^\{2\}values\. These failures are largely attributable to PyVision’s implementation strategy, which is primarily optimized around MAP\-based fitting and does not consistently translate to stable MCMC\-based inference\.

## Appendix EExpert Tools in Detail

### E\.1Distribution Fitting

1. 1\.CalculateMoments: Computes the mean, variance, skewness, and excess kurtosis of the input data\. Returns a JSON artifact with a plain\-language interpretation to guide distribution selection, including symmetry hints \(e\.g\., right\-skewed data suggests Gamma, Lognormal, or Weibull families\) and tail\-weight hints \(e\.g\., leptokurtic data suggests Student\-tt, Cauchy, or Laplace families\)\.
2. 2\.Histogram: Plots a histogram of the empirical data with the fitted distribution’s probability density function \(PDF\) overlaid\. Handles both single distributions and mixtures by summing component PDFs weighted by their mixture weights\. Provides an immediate visual check of whether the model captures the overall shape, modality, and spread of the data\. When no model has been fitted yet, produces a raw histogram without an overlay\.
3. 3\.SegmentDistributionsAndCalculateMoments: Segments the data into a specified number of mixture components using a Gaussian Mixture Model \(GMM\), then computes per\-component moments \(mean, variance, skewness, kurtosis\)\. Produces both a segmentation image with a total mixture overlay and a JSON summary of per\-component statistics with distribution family hints\. Most useful when the histogram appears multimodal\.
4. 4\.QQPlot: Generates a Quantile\-Quantile \(Q\-Q\) plot comparing empirical data quantiles to theoretical quantiles from the currently fitted distribution\. Linearity indicates a good fit; S\-shaped curvature signals tail mismatch; one\-sided curvature suggests skew; and sharp tail departures may indicate outliers or heavier tails than the model captures\.
5. 5\.PlotTailsTransform: Produces log\-log and semi\-log complementary CDF \(CCDF\) plots to diagnose tail behavior\. A straight line on the log\-log plot indicates power\-law or Pareto\-type heavy tails, while a straight line on the semi\-log plot indicates exponential decay\. Useful for distinguishing heavy\-tailed from light\-tailed distributions when the histogram alone is ambiguous\.
6. 6\.ProbabilityPlot: Generates a probability plot comparing the empirical CDF to the fitted distribution’s theoretical CDF\. A consistent horizontal shift indicates a mis\-specified location parameter; a slope mismatch indicates a scale misfit; and systematic tail deviations suggest distributional misfit\. Also reports a Kolmogorov\-Smirnov \(KS\) statistic for quantitative goodness\-of\-fit assessment\.

### E\.2Time Series

Table 9:Per\-sample tool coverage: percentage of samples in which each Expert tool is invoked at least once during the iterative refinement loop\.1. 1\.GetDominantPeriod: Extracts the dominant period from the time series using Fast Fourier Transform \(FFT\) analysis\. Returns a plain\-text summary of the detected period\. Most useful when Periodic or PeriodicComplex kernels are under consideration and the period has not yet been numerically determined\. The result is available in the subsequent feedback iteration\.
2. 2\.FitVsActuals: Produces a visual overlay of the Gaussian Process \(GP\) fit on the raw time series data\. Essential for visually assessing whether the model adequately captures the underlying trend and seasonality while appropriately discounting noise\. Falls back to a raw series plot if no model has been fitted yet\.
3. 3\.FitVsActualsWithResidualsDistribution: Generates a combined plot showing the GP fit overlaid on the observed time series alongside the distribution of residuals\. Used to assess whether residuals resemble white noise; a broadly normal residual distribution is indicative of a well\-specified model\. Falls back to a raw series plot if no model has been fitted yet\.
4. 4\.ResidualsAutoCorrelationPlot: Produces an Autocorrelation Function \(ACF\) plot of the model residuals to check for temporal independence\. Significant spikes above the confidence band indicate that the model is failing to capture some latent structure in the data\. Falls back to a raw series plot if no model has been fitted yet\.
5. 5\.ResidualsAutoCorrelationScore: Applies the Ljung\-Box statistical test to assess residual independence numerically\. App\-value greater than0\.050\.05indicates that the residuals are consistent with white noise, suggesting an adequate model fit\. Returns a plain\-text summary and requires a fitted model to produce meaningful results\.

#### E\.2\.1Feedback Modality and Adaptive Selection

The static and dynamic toolkits differ fundamentally in their feedback modality\. The expert toolkit provides narrow, targeted signals—such as returning a single scalar for dominant frequency—which allows the VLM to rapidly integrate findings and converge quickly on structurally simple datasets \(Easy, Chirp\)\.

In contrast, the dynamic toolkit generates comprehensive, multi\-panel visualizations that synthesize frequency estimation, residual analysis, and harmonic structure into a single artifact\. While this rich feedback demands more visual interpretation from the VLM, it is critical for diagnosing complex signals where isolated scalars are insufficient\.

As shown in Table[9](https://arxiv.org/html/2606.00384#A5.T9), the VLM dynamically adapts its tool selection to the underlying problem complexity\. Simple goodness\-of\-fit overlays dominate the early steps for clean periodic datasets\. However, when applied to the Hard dataset, the agent heavily prioritizes autocorrelation and residual diagnostics to detect more subtle, persistent temporal dependencies\.

#### E\.2\.2Comparison with PyVision Baseline

Apart of choice of modalities, and complexity of serving the tool \(via panels\), both dynamic\-toolkit and expert\-toolkit lean into visual feedback\. On the other hand, PyVision also generates dynamic tools on the fly, but they diverge significantly in diagnostic focus, output density, and tool lifecycle\.

##### Diagnostic Category Coverage\.

As detailed in Table[10](https://arxiv.org/html/2606.00384#A5.T10), both systems heavily utilize Lag/Autocorrelation and Frequency Estimation \(∼\{\\sim\}85–94%\)\. However, VESTA’s dynamic toolkit exhibits a massive preference for Residual Analysis \(91\.9% vs\. PyVision’s 11\.9%\) and exclusively employs Statistical Tests \(28\.6% vs\. 0%\)\. Because residual structure is the primary indicator of GP kernel misspecification, this gap directly impacts model selection quality\.

Table 10:Diagnostic category coverage \(% of samples with≥\\geq1 tool in category\)\. PV = PyVision\.

## Appendix FAnalysis of Generated Tools Not in the Expert Toolkit

### F\.1Distribution Fitting

1. 1\.diagnostic\_fit\_checks:Naming a concrete model family \(gaussian, gamma, lognormal, Pareto, Weibull, etc\.\) and trying it on the data\. Each family encodes different assumptions\. These tools allow for typically allow for a visual comparison of multiple model families at once\. Occasionally, we observe some single use model fitting\.
2. 2\.information\_criteria:Numerical scores that rank competing fits while penalizing model complexity\. Beyond simply fitting and visualizing models, AIC and BIC provide quantitative fit metrics\. Table 11:Analysis of functions inVesta\-generated tools that are not contained in the expert toolkit for Distribution Fitting\.
3. 3\.mle\_fitting:Maximum likelihood estimation: choosing the parameter values \(μ\\mu,σ\\sigma, shape, scale, …\) that maximize the probability of observing the data under the chosen family\. This is the how you actually of fit models, distinct from what models we want to test indiagnostic\_fit\_checks\. MLE gives you the canonical “best” parameters under a given family, plus a likelihood value that feeds directly into AIC/BIC\.
4. 4\.mean\_excess\_plotPlots the conditional expectation𝔼​\[X−u​∣X\>​u\]\\mathbb\{E\}\[X\-u\\mid X\>u\]against thresholduu\. For the Generalized Pareto distribution this function is linear inuu, so a straight line in the upper tail signals a GPD\-like tail and tells you where the “extreme regime” begins\. This is a tail\-diagnostic that complementsdiagnostic\_fit\_checks\. These test help itVestadecide if the tail behaves as that family predicts and at what threshold the heavy\-tail regime kicks in\. The 36% rate on Astro vs\.∼\\sim5% elsewhere is a sign the model recognizes astronomical data as tail\-dominated\.
5. 5\.hill\_estimatorEstimates the tail indexα\\alphaof a heavy\-tailed distribution from the largestkkorder statistics, giving a concrete number for “how heavy” the tail is\. A Hill plot \(α^\\hat\{\\alpha\}vs\.kk\) lets you check stability and pick a sensible threshold\. This refines a Pareto/power\-law fit by pinning down its single most important parameter, and serves as a sanity check on MLE\-fit tail parameters — if MLE on a full\-data Pareto fit and a Hill estimate on the upper tail disagree, the body of the distribution is contaminating your tail estimate\.
6. 6\.shapiro\_wilkA formal hypothesis test for whether data come from a normal distribution\. If Shapiro\-Wilk rejects normality strongly, that rules out the normal family indiagnostic\_fit\_checks\.
7. 7\.box\_coxA parametric family of power transformsy=\(xλ−1\)/λy=\(x^\{\\lambda\}\-1\)/\\lambdathat searches for theλ\\lambdamaking the transformed data closest to normal and can be useful when working with exotic, heavy\-tailed distributions\.

Table 12:Analysis of functions inVesta\-generated tools that are not contained in the expert toolkit for Distribution Fitting\.
### F\.2Time Series

1. 1\.density\_visualization:Overlays a histogram with a kernel density estimate \(KDE\) to give a non\-parametric picture of the marginal distribution of a time series\. This is typically the first diagnosticVestareaches for: before committing to any parametric family, the histogram/KDE pair reveals modality, skew, and gross tail behavior at a glance\. The KDE bandwidth controls the bias–variance tradeoff in smoothing; narrower bandwidths reveal fine structure while wider ones expose the global shape\.
2. 2\.distribution\_shape\_summary:Computes the first four standardized moments \(mean, variance, skewness, and excess kurtosis\) along with the coefficient of variationCV=σ/μ\\text\{CV\}=\\sigma/\\mu\. These scalars compress the shape of an entire distribution into a small, comparable summary\. Skewness distinguishes asymmetric processes \(e\.g\. log\-normal dwell times\) from symmetric ones\. Kurtosis flags heavy tails before any formal test is run\. This is the numerical complement todensity\_visualization\.
3. 3\.moving\_average:Computes a rolling meanx¯t=1w​∑i=0w−1xt−i\\bar\{x\}\_\{t\}=\\frac\{1\}\{w\}\\sum\_\{i=0\}^\{w\-1\}x\_\{t\-i\}over a sliding window of widthww\. This is low\-pass filter suppresses high\-frequency noise and reveals the underlying trend or long\-timescale signal\.Vestauses it both as a standalone smoother and as a preprocessing step before stationarity tests or ARIMA fitting, since many parametric models assume a trend\-free series\.
4. 4\.regime\_switching\_analysis:A grouped capability covering Gaussian Mixture Models \(GMMs\), bimodality and multimodality detection,kk\-means clustering on temporal states, Hidden Markov Models \(HMMs\), and log\-normal fits to dwell times\. Together these tools test the amount of time a signal spends in qualitatively distinct states?, GMMs andkk\-means partition the amplitude space into clusters; HMMs additionally impose a Markov transition structure on state sequences; dwell\-time fits characterize how long the process stays in each regime\. This matters for telegraph\-noise or switching processes where a single stationary distribution is a category error\. The high rate onHardreflects that harder series are more likely to be genuine multi\-state processes\.
5. 5\.heavy\_tail\_diagnostics:A grouped capability covering Pareto/power\-law identification, Generalized Extreme Value \(GEV\) fits, Hill estimator calls, and mean\-excess plots; see the analogous entries in the distribution\-fitting appendix for mathematical detail on the latter two\. In the time series context,Vestaapplies these when the marginal distribution of the series itself \(rather than residuals\) appears heavy\-tailed, e\.g\. burst\-like signals or processes with rare but extreme excursions\.
6. 6\.qq\_plot:Plots empirical quantilesFn−1​\(p\)F\_\{n\}^\{\-1\}\(p\)against theoretical quantilesFθ−1​\(p\)F\_\{\\theta\}^\{\-1\}\(p\)for a candidate familyFθF\_\{\\theta\}\. These are typically used to measure the distribution of residuals\.
7. 7\.shapiro\_wilk:A formal hypothesis test of the null that the data are drawn from a normal distribution, with test statisticW=\(∑ai​x\(i\)\)2/∑\(xi−x¯\)2W=\(\\sum a\_\{i\}x\_\{\(i\)\}\)^\{2\}/\\sum\(x\_\{i\}\-\\bar\{x\}\)^\{2\}wherex\(i\)x\_\{\(i\)\}are order statistics andaia\_\{i\}are weights derived from the expected normal order statistics\. In the time series context,Vestatypically applies this to residuals after trend/seasonality removal to verify that the generating process is Gaussian, which is a prerequisite for classical ARIMA inference\.
8. 8\.stationarity\_tests:A grouped capability covering the Augmented Dickey–Fuller \(ADF\) test for a unit root and the KPSS test for level/trend stationarity\. These are complementary\. ADF tests the null of a unit root while KPSS tests the null of stationarity, so using both guards against the failure mode of each test individually\. Stationarity is a prerequisite for most parametric time series models; this capability gates whetherVestaproceeds to ARIMA fitting or first applies differencing or detrending\.
9. 9\.differencing:Computes the lag\-dddifference∇dxt=xt−xt−d\\nabla^\{d\}x\_\{t\}=x\_\{t\}\-x\_\{t\-d\}, most commonlyd=1d=1, the first difference\. This is the standard cure for a unit root: if the ADF test fails to reject non\-stationarity, first\-differencing typically removes a stochastic trend and produces a stationary series suitable for ARIMA fitting\. Seasonal differencing \(d=sd=s, the period\) removes periodic non\-stationarity\.Vestauses this as a preprocessing step that is directly triggered by stationarity test results\.
10. 10\.arima:Fits an Autoregressive Integrated Moving Average modelϕ​\(B\)​\(1−B\)d​xt=θ​\(B\)​ϵt\\phi\(B\)\(1\-B\)^\{d\}x\_\{t\}=\\theta\(B\)\\epsilon\_\{t\}, whereBBis the backshift operator,ppautoregressive terms capture dependence on recent values,dddifferences remove trends, andqqmoving\-average terms model the residual autocorrelation structure\. ARIMA is the workhorse parametric model for univariate stationary time series; it is typically applied after stationarity tests confirm \(or differencing achieves\) stationarity, and model order\(p,d,q\)\(p,d,q\)is selected via AIC/BIC\.
11. 11\.change\_point\_detection:Identifies timesτ\\tauat which the statistical properties of the series, mean, variance, or autocorrelation structure, shift abruptly\. Unlikeregime\_switching\_analysis, which models soft probabilistic membership in states, change\-point methods partition the timeline into contiguous segments and locate the breakpoints\. This is appropriate when the series undergoes irreversible structural shifts \(e\.g\. a one\-time intervention\) rather than recurrent toggling between states\.
12. 12\.information\_criteria:Computes AIC=2​k−2​ℓ^=2k\-2\\hat\{\\ell\}and BIC=k​ln⁡n−2​ℓ^=k\\ln n\-2\\hat\{\\ell\}, wherekkis the number of free parameters,nnthe sample size, andℓ^\\hat\{\\ell\}the maximized log\-likelihood\. These penalize model complexity to prevent overfitting and give you a principled basis for comparing models of different orders \(e\.g\. ARIMA\(1,1,1\)\(1,1,1\)vs\. ARIMA\(2,1,2\)\(2,1,2\)\) or different parametric families on the same data\. AIC favors predictive accuracy while BIC applies a stronger complexity penalty and is consistent for model selection\.
13. 13\.mle\_fitting:See the analogous entry in the distribution\-fitting appendix\. In the time series context, MLE is used both for fitting marginal distributions \(e\.g\. the amplitude distribution of a stationary series\) and for estimating ARIMA parameters via conditional or exact likelihood\. The resulting likelihood value feeds directly into AIC/BIC for model comparison\.
14. 14\.kolmogorov\_smirnov:A nonparametric goodness\-of\-fit test based on the supremum distanceDn=supx\|Fn​\(x\)−F0​\(x\)\|D\_\{n\}=\\sup\_\{x\}\|F\_\{n\}\(x\)\-F\_\{0\}\(x\)\|between the empirical CDF and a hypothesized distributionF0F\_\{0\}\. Unlike Shapiro–Wilk, the KS test is distribution\-free and applies to any fully specifiedF0F\_\{0\}, making it useful whenVestawants to test an arbitrary fitted family rather than only normality\. It is less powerful than Anderson–Darling in the tails but gives a straightforward visual diagnostic via the CDF overlay\.
15. 15\.log\_transform:Appliesyt=log⁡\(xt\)y\_\{t\}=\\log\(x\_\{t\}\)to stabilize variance in series whose standard deviation grows with the mean, a signature of log\-normal or multiplicative processes\. Log\-transforming a series with exponential growth or multiplicative noise often produces an approximately Gaussian, stationary series amenable to ARIMA\. It is also a prerequisite for power\-law tail analysis, where log–log linearity of the CCDF is the diagnostic signature\.
16. 16\.jarque\_bera:A moment\-based normality test with statisticJ​B=n6​\(γ12\+γ224\)JB=\\frac\{n\}\{6\}\\left\(\\gamma\_\{1\}^\{2\}\+\\frac\{\\gamma\_\{2\}^\{2\}\}\{4\}\\right\), which is asymptoticallyχ2​\(2\)\\chi^\{2\}\(2\)under normality\. Because it is constructed directly from skewness and kurtosis, it is more powerful than the KS test against symmetric heavy\-tailed alternatives \(largeγ2\\gamma\_\{2\}\) and skewed alternatives \(large\|γ1\|\|\\gamma\_\{1\}\|\)\. In practiceVestauses it as a fast scalar complement to the fulldistribution\_shape\_summary\.
17. 17\.ecdf:Plots the empirical cumulative distribution functionFn​\(x\)=1n​∑i=1n𝟏​\[xi≤x\]F\_\{n\}\(x\)=\\frac\{1\}\{n\}\\sum\_\{i=1\}^\{n\}\\mathbf\{1\}\[x\_\{i\}\\leq x\]\. The ECDF is a non\-parametric, unsmoothed estimate of the true CDF and serves as the ground truth against which parametric fits are visually compared\. It is the natural complement to the KS test: the test statistic is the maximum vertical gap between the ECDF and the hypothesized CDF curve\.

## Appendix GRuntime ofVesta

Table 13:Runtimes by Domain, Dataset, Mode, and Model
## Appendix HPrompts

This section contains all the prompts used by our method and baselines\.

VestaTime Series Critique PromptYou are the VLM critic ofVesta\. Inspect the latest fit of the current Gaussian\-process PyMC model and decide whether to refine the kernel composition further or to terminate the iterative refinement loop\. The prompt covers stationary periodic structure as well as non\-stationary \(warped periodic\) wave signals\.Goal:Determine whether the current GP kernel composition correctly captures the trend, seasonality, and non\-stationarity of the series, and whether the priors are well\-scaled for the observed amplitude and period\.Current Model Information•Current PyMC code:\{current\_pymc\_code\}•Current kernels:\{current\_kernels\}•Previously tested kernel combinations:\{tested\_kernels\}•History of previously executed tools:\{tool\_history\}Iteration Summaries:\{summaries\}Important Constraints•Evaluate fit quality based on the most recent iteration’s tool output and visual feedback\.•Take the trajectory of past iterations and their metric values into account\.•You may revisit a previously tested kernel combination, but you must use sufficiently different priors\.•Do not invoke the same tool more than twice across the run\.•Combine kernels by addition only \(cov = cov1 \+ cov2\)\.•warped\_periodicsubsumesperiodic; do not combine them\.Stationary vs\. Non\-Stationary Wave Check:\[ visual rule for deciding whether cycle spacing is constant \(stationary\) or sweeping \(non\-stationary\), and what kernel that implies \]Kernel Family Guide:\[ valid kernel names, when each kernel applies, and how kernels may be combined are listed here \]Prior Specification Rules:\[ recommended prior families and scales for length\-scale, amplitude, period, warp coefficients, and observation noise are listed here, with anti\-patterns for misuse ofpm\.Deterministicand undefined initializer variables \]Proposal Diversity Rules:\[ the primary proposal must be the best guess; the remaining proposals must explore structurally different kernel combinations, including at least one with and one withoutwarped\_periodicwhen non\-stationary periodicity is suspected \]Diagnostic Results:\{tool\_outputs\}Task:If the current fit is satisfactory, terminate the loop early by including the literal tokenCOMPLETEin your description\. Otherwise, propose exactly\{num\_proposals\}diverse revised models following the rules above\. State explicitly whether the stationary\-vs\-non\-stationary wave check passes and, if non\-stationary, which direction the frequency sweeps\.Output Format:Returnonlya valid JSON object withdescriptionandproposals\. Each proposal must includekernelsandpriors\. Never output anything outside the JSON object\.Figure 7:Critique\-stage prompt used byVestafor time series modeling\.VestaDistribution Fitting Critique PromptYou are the VLM critic ofVesta\. Inspect the latest fit of the current PyMC model and decide whether to refine it further or to terminate the iterative refinement loop\.Goal:Determine whether the current distribution model correctly captures the shape, location, scale, and tail behavior of the data\.Current Model Information•Current PyMC code:\{current\_pymc\_code\}•Current distribution family:\{current\_family\}•Previously tested distribution families:\{tested\_families\}•History of previously executed tools:\{tool\_history\}Iteration Summaries:\{summaries\}Important Constraints•Evaluate fit quality based on the most recent iteration’s tool output and visual feedback\.•Take the trajectory of past iterations and their metric values into account\.•You may revisit a previously tested family, but you must use sufficiently different priors\.•Do not invoke the same tool more than twice across the run\.•Prefer interpretable distributions and avoid unnecessary complexity\.Model Fit Diagnostics:\[ shape capture, peak alignment, tail behavior, mixture detection, and overfitting checks are listed here \]Distribution Family Guide:\[ valid single families, mixture syntax, parametrization reference for each family, support\-validation rules, and family\-specific prior constraints are listed here \]Prior Specification Rules:\[ instructions for reading the histogram to extract peak/spread/support/tails, translating those observations into hyperparameter values, and the variable naming convention for single and mixture priors are listed here \]Proposal Diversity Rules:\[ the primary proposal must be the best guess; the remaining proposals must be structurally different alternatives covering both single distributions and mixtures \]Diagnostic Results:\{tool\_outputs\}Task:If the current fit is satisfactory, terminate the loop early by including the literal tokenCOMPLETEin your description\. Otherwise, propose exactly\{num\_proposals\}diverse revised models following the rules above\.Output Format:Returnonlya valid JSON object withdescriptionandproposals\. Each proposal must includedistribution\_family,is\_mixture, andpriors\. Never output anything outside the JSON object\.Figure 8:Critique\-stage prompt used byVestafor distribution fitting\.VestaTool Selection PromptYou are evaluating the current PyMC model fit\. The latest fit is shown in the attached plot\.Tool Registry:\{tool\_descriptions\}Decision Rules•Select exactly one tool from the registry that best diagnoses the current fit\.•Prefer tools that probe the suspected weakness of the current model:\{model\_weakness\_examples\}\.•If no existing tool in the registry is adequate, return thegenerate\_new\_toolaction and describe what the new tool should compute; the tool\-creation step will then code it\.•You may not invoke the same tool more than twice across the run\.Output Format:Returnonlya valid JSON object naming the selected tool and any arguments it requires\.Domain\-specific values for the placeholders above:•Distribution fitting\-\{model\_weakness\_examples\}=peak misalignment, tail mismatch, multimodality, support violation\.•Time series modeling\-\{model\_weakness\_examples\}=missed period, residual structure, amplitude envelope, frequency sweep, anomalies\.Figure 9:Tool\-selection prompt used by the Generate\-Tools stage ofVesta\.VestaTool Creation PromptYou are creating a new diagnostic tool that will be added to the tool registry ofVestaand may be reused in later iterations of this run\. The tool is a Python function executed in a sandbox with access to:numpy as np,pandas as pd,scipy\.stats,matplotlib\.pyplot as plt,sklearn\{extra\_imports\}\.Function Signature•Define a single function:def tool\(data: \{data\_handle\}, map\_estimate: dict, family\_or\_kernel: str\) \-\> str•The return value is the path to a saved figure \(PNG\) or a short textual summary that the VLM critic can read\.Implementation Rules•Begin directly with thedef tool\(\.\.\.\)block\. Donotinclude import statements; required modules are already in scope\.•Usemap\_estimateto construct a comparison between the observed data and the current model\.•Save any figure withplt\.savefig\(\.\.\.\)and return its path; do not callplt\.show\(\)\.•The tool must run end\-to\-end without user interaction\.•\{domain\_specific\_code\_rule\}Output Format:Returnonlya JSON\-parseable dict with key"code"containing the function as a single\-line string with\\nfor newlines\.Domain\-specific values for the placeholders above:•Distribution fitting\-\{extra\_imports\}= \(none beyond the default sandbox\);\{data\_handle\}=np\.ndarray;\{domain\_specific\_code\_rule\}= “Visualize the empirical and model\-implied densities or transformed\-data summaries; donotfit a new PyMC model inside the tool\.”•Time series modeling\-\{extra\_imports\}=, statsmodels\.api as sm, scipy\.signal;\{data\_handle\}=pd\.Series with a DatetimeIndex;\{domain\_specific\_code\_rule\}= “Use the time index for any plot’s x\-axis and align predictions withdata\.index; donotre\-fit the GP inside the tool\.”Figure 10:Tool\-creation prompt used by the Generate\-Tools stage ofVesta\.VestaSummarization PromptYou are summarizing the latest iteration of a\{run\_kind\}run so that the next iteration ofVestacan reason over the full refinement trajectory without consuming the raw intermediate outputs\.Inputs Available•The visual description produced for the current iteration\.•The best PyMC model code and its fitted parameters \(map\_estimate\)\.•The metric value and any auxiliary fit statistics\.•The output \(image or text\) produced by the selected or newly created tool\.What the Summary Must Contain•A short textual description of what the data look like and what structure the current best model captures\.•The chosen\{model\_field\_name\}and its key fitted parameters as concise key–value pairs\.•The metric value for the current best model\.•One or two sentences describing what the tool output revealed, focusing on remaining mismatch between the model and the data\.•\{description\_specialization\}Style Rules•Be concrete: name parameters and numeric values rather than describing them abstractly\.•Donotre\-derive the model or restate prompt instructions\.•Keep the summary self\-contained; it replaces the raw per\-iteration outputs in subsequent iterations\.Output Format:Return a JSON object with fieldsdescription,\{model\_field\_name\},map\_estimate,metric, andtool\_output\_summary\.Domain\-specific values for the placeholders above:•Distribution fitting\-\{run\_kind\}=distribution fitting;\{model\_field\_name\}=distribution\_family;\{description\_specialization\}= “Mention any tail, peak, or modality features the current model fails to reproduce\.”•Time series modeling\-\{run\_kind\}=time series modeling;\{model\_field\_name\}=kernels;\{description\_specialization\}= “Mention any trend, seasonality, or non\-stationarity features the current kernels fail to reproduce\.”Figure 11:Prompt used by the Summarize stage ofVesta\.BoxLM Proposal Prompt for Time Series ModelingYou are a brilliant statistician modeling a time series dataset using Gaussian processes\. Your job is to come up with a GP model that explains the time series by writing a PyMC probabilistic program\. In the first round you receive only the dataset description; in later rounds you also receive the previous program along with the critic’s hypotheses and synthesis\.Inputs•Dataset description:\{dataset\_description\}•Column description:\{column\_description\}•Previous round’s hypotheses:\{prev\_hypotheses\}•Previous round’s synthesis:\{prev\_synthesis\}Kernel Family Guide•pm\.gp\.cov\.Linear\- persistent upward or downward drift\.•pm\.gp\.cov\.Periodic\- smooth repeating sinusoidal patterns\.•pm\.gp\.cov\.ExpQuad\(RBF\) \- extremely smooth curves\.•pm\.gp\.cov\.Matern52\- rough or jagged paths\.•Combine kernels by addition for independent structures, by multiplication for modulation\. Do not combine RBF and Matern in the same model\.PyMC API Rules•Importpymc\(notpymc3\)\. Do not usesdas an argument; usesigmainstead\.•Use the variable namey\_obsfor the likelihood\.•Index columns explicitly fromobserved\_data; never useobserved\_data\.indexto grab raw data\.•The model must use Gaussian processes \(pm\.gp\); plain regression is not allowed\.•Do not usepm\.Constantfor parameters and do not passdims="obs\_id"for the GP marginal likelihood\.Required Function Skeleton•Define a functiondef gen\_model\(observed\_data\)that loads data viaobserved\_data\["column\_name"\]\.values\(never a placeholder array literal\)\.•Before the model block, normalize the time axis to\[0, 1\]and center the target values\.•Insidewith pm\.Model\(\) as model:declare the kernel priors, build the covariance, definegp = pm\.gp\.Marginal\(cov\_func=cov\), and set the likelihoody\_obs = gp\.marginal\_likelihood\(’y\_obs’, \.\.\., sigma=sigma\)\.•Sample withtrace = pm\.sample\(1000, tune=500, target\_accept=0\.90, chains=3, cores=1, random\_seed=rng1, idata\_kwargs=\{"log\_likelihood": True\}\)\.•Then runposterior\_predictive = pm\.sample\_posterior\_predictive\(trace, random\_seed=rng2, return\_inferencedata=False\)\.•Return\(model, posterior\_predictive, trace\)in that order\.Response Structure•First, sketch a high\-level GP model for the data and explicitly cite which previous\-round hypotheses \(if any\) the sketch addresses\.•Then write the full program in a Python code block, with comments on lines that address specific hypotheses\.Figure 12:Proposal prompt used by the BoxLM baseline for time series modeling\.BoxLM Proposal Prompt for Distribution FittingYou are a brilliant statistician modeling a dataset\. Your job is to come up with a generative model that explains the true data by writing a PyMC probabilistic program\. In the first round you receive only the dataset description; in later rounds you also receive the previous program along with the critic’s hypotheses and synthesis\.Inputs•Dataset description:\{dataset\_description\}•Column description:\{column\_description\}•Previous round’s hypotheses:\{prev\_hypotheses\}•Previous round’s synthesis:\{prev\_synthesis\}PyMC API Rules•Importpymc\(notpymc3\)\. Do not usesdas an argument; usesigmainstead\.•Use the variable namey\_obsfor the likelihood\.•Index columns explicitly fromobserved\_data; never useobserved\_data\.indexto grab raw data\.•Wrap every input feature in apm\.MutableDatacontainer withdims="obs\_id", and passdims="obs\_id"to the likelihood as well\.•Do not usepm\.Constantfor parameters and do not use Gaussian processes for this task\.•For Student\-T, givenua continuous positive prior such aspm\.Exponentialorpm\.HalfNormal\.Required Function Skeleton•Define a functiondef gen\_model\(observed\_data\)that loads data viaobserved\_data\["column\_name"\]\.values\(never a placeholder array literal\)\.•Insidewith pm\.Model\(\) as model:declare theMutableDatacontainers, the priors, and the likelihoody\_obs\.•Sample withtrace = pm\.sample\(1000, tune=500, target\_accept=0\.90, chains=3, cores=1, random\_seed=rng1, idata\_kwargs=\{"log\_likelihood": True\}\)\.•Then runposterior\_predictive = pm\.sample\_posterior\_predictive\(trace, random\_seed=rng2, return\_inferencedata=False\)\.•Return\(model, posterior\_predictive, trace\)in that order\.Response Structure•First, sketch a high\-level probabilistic program for the data and explicitly cite which previous\-round hypotheses \(if any\) the sketch addresses\.•Then write the full program in a Python code block, with comments on lines that address specific hypotheses\.Figure 13:Proposal prompt used by the BoxLM baseline for distribution fitting\.BoxLM Critic PromptYou are a brilliant statistician specializing in critiquing and proposing revisions of probabilistic models\. Your colleague has produced a list of PyMC programs that hypothesize\{model\_class\}for the data; each program has been fit, plotted against the actual data, and scored with ELPD\-LOO \(higher is better\)\. Your job is to figure out whether the proposed programs are consistent with the actual data by comparing the posterior predictive plots against the actual data, and to provide modeling\-related guidance for the next round of revision\.Inputs•Dataset description:\{dataset\_description\}•Column description:\{column\_description\}•Expert context:\{expert\_context\}•Programs from this round, each with its ELPD\-LOO score and a posterior\-predictive plot overlaid on the actual data\.•Hypotheses and synthesis from previous rounds \(if any\), which may be revised, removed, or extended\.Instructions•Look at every posterior\-predictive plot and explicitly describe the visual properties of the model versus the data\.•State natural\-language hypotheses for any discrepancies between the proposed programs and the actual data, prioritizing the larger discrepancies\.•Using the programs and their ELPD\-LOO scores, write a synthesis of which strategies and modeling choices did and did not work\. The synthesis is passed back to your colleague, so make it as informative as possible\. Restrict suggestions to modeling; do not suggest collecting more data\.•If hypotheses or a synthesis from previous rounds are provided, revise them based on the new evidence and briefly explain any additions or removals\.Output Format
Respond with two fenced blocks in this exact order:```
‘‘‘ Hypotheses
    * Hypothesis
    * Hypothesis
‘‘‘

‘‘‘Synthesis
    (Your synthesis here)
‘‘‘
```

Domain\-specific values•For distribution fitting,\{model\_class\}is*generative models*\.•For time series modeling,\{model\_class\}is*Gaussian process models*\.Figure 14:Critic prompt used by the BoxLM baseline for both distribution fitting and time series modeling\.PyVision Agent PromptYou are an agent; please keep going until the user’s query is completely resolved before ending your turn and yielding back to the user\. Only terminate your turn when you are sure that the problem is solved\. Solve the following problem step by step\. You now have the ability to selectively write executable Python code to enhance your reasoning process\. The Python code will be executed by an external sandbox\.You MUST plan extensively before each function call, and reflect extensively on the outcomes of the previous function calls\. Do not do this entire process by making function calls only, as this can impair your ability to solve the problem and think insightfully\.For all the provided images, in order, the*i*\-th image has already been read into the global variableimage\_clue\_iusingPIL\.Image\.open\(\)\. When writing Python code, you can directly use these variables without needing to read them again\.Since you are dealing with the vision\-related task, you MUST use the Python tool \(e\.g\., thematplotliblibrary\) to analyze or transform the image whenever it could improve your understanding or aid your reasoning\. This includes but is not limited to zooming in, rotating, adjusting contrast, computing statistics, or isolating features\. When you usematplotlibto visualize data or further process images, useplt\.show\(\)to display them; there is no need to save them\. Do not use image\-processing libraries likecv2orPIL\. If you want to check the value of a variable, you MUST useprint\(\)to check it\.Code\-execution Format
The output \(wrapped in<interpreter\>output\_str</interpreter\>\) is returned to aid your reasoning\. Each Python snippet must be a complete script including necessary imports, wrapped exactly as:```
<code>
‘‘‘python
code snippet
‘‘‘
</code>
```

Image Information
Image Width: \{width\}; Image Height: \{height\}User Query
\{query\}Final Answer Format
The last part of your response must be in the following format:```
<answer>
\boxed{’The final answer goes here.’}
</answer>
```

Domain\-specific values•For distribution fitting,\{query\}instructs the agent to propose and fit a PyMC probabilistic program for the empirical distribution shown in the histogram image, using the same PyMC coding guidelines as theVestaagent\.•For time series modeling,\{query\}instructs the agent to propose and fit a PyMC Gaussian process program for the series shown in the time series plot, using the same PyMC coding guidelines as theVestaagent\.Figure 15:Agent prompt used by the PyVision baseline for both distribution fitting and time series modeling\.
## Appendix ILimitations

WhileVestademonstrates strong performance acrossDawn, several limitations remain\. First,Vestarelies on VLMs to interpret multi\-panel diagnostic visualizations, and current VLMs struggle to fully unpack and reason over complex visual outputs\. This bottleneck directly limits the gains from dynamic tool creation, as the most sophisticated generated tools contain multi\-panel images that compose multiple functions\. This gap likely accounts for the remaining performance differential betweenVesta\-Dynamic andVesta\-Expert\.

Second, the iterative nature ofVestaincurs non\-trivial computational cost, as each iteration requires fitting up toppPyMC models via MAP estimation\. This limits scalability to very large datasets\.

Finally, whileDawncovers two domains central to scientific modeling and includes real\-worldAstrotasks, the benchmark is synthetically generated\. Performance on fully observational scientific datasets, where ground\-truth parameters are unknown, remains an important avenue for future evaluation\.