WO2024063907A1

WO2024063907A1 - Modelling causation in machine learning

Info

Publication number: WO2024063907A1
Application number: PCT/US2023/031000
Authority: WO
Inventors: Wenbo GONG; Cheng Zhang; Nick Pawlowski; Joel Jennings; Karen FASSIO; Marife DEFANTE; Steve Thomas; Alice HORAN; Chao Ma; Matthew ASHMAN; Agrin HILMKIL
Original assignee: Microsoft Technology Licensing LLC
Current assignee: Microsoft Technology Licensing LLC
Priority date: 2022-09-21
Filing date: 2023-08-24
Publication date: 2024-03-28
Anticipated expiration: 2025-03-21

Abstract

A method comprising: sampling a temporal causal graph from a temporal graph distribution specifying probabilities of directed causal edges between different variables of a feature vector at a present time step, and from one variable at a preceding time step to another variables at the present time step. Based on this there are identified: a present parent which is a cause of the selected variable in the present time step, and a preceding parent which is a cause of the selected variable from the preceding time step. The method then comprises: inputting a value of each identified present and preceding parent into a respective encoder, resulting in a respective embedding of each of the present and preceding parents; combining the embeddings of the present and preceding parents, resulting in a combined embedding; inputting the combined embedding into a decoder, resulting in a reconstructed value of the selected variable.

Description

MODELLING CAUSATION IN MACHINE LEARNING BACKGROUND Neural networks are used in the field of machine learning and artificial intelligence (AI). A neural network comprises plurality of nodes which are interconnected by links, sometimes referred to as edges. The input edges of one or more nodes form the input of the network as a whole, and the output edges of one or more other nodes form the output of the network as a whole, whilst the output edges of various nodes within the network form the input edges to other nodes. Each node represents a function of its input edge(s) weighted by a respective weight, the result being output(s) on its output edge(s). The weights can be gradually tuned based on a set of training data so as to tend towards a state where the output of the network will output a desired value for a given input.Typically the nodes are arranged into layers with at least an input and an output layer. A “deep” neural network comprises one or more intermediate or “hidden” layers in between the input layer and the output layer. The neural network can take input data and propagate the input data through the layers of the network to generate output data. Certain nodes within the network perform operations on the data, and the result of those operations is passed to other nodes, and so on. Each node is configured to generate an output by carrying out a function on the values input to that node. The inputs to one or more nodes form the input of the neural network, the outputs of some nodes form the inputs to other nodes, and the outputs of one or more nodes form the output of the network. At some or all of the nodes of the network, the input to that node is weighted by a respective weight. A weight may define the connectivity between a node in a given layer and the nodes in the next layer of the neural network. A weight can take the form of a scalar or a probabilistic distribution. When the weights are defined by a distribution, as in a Bayesian model, the neural network can be fully probabilistic and captures the concept of uncertainty. The values of the connections between nodes may also be modelled as distributions. The distributions may be represented in the form of a set of samples or a set of parameters parameterizing the distribution (e.g. the mean μ and standard deviation σ or variance σ²).The network learns by operating on data input at the input layer, and adjusting the weights applied by some or all of the nodes based on the input data. There are different learning approaches, but in general there is a forward propagation through the network from input layer to output layer, a calculation of an overall error, and a backward propagation of the error through the network from output layer to input layer. In the next cycle, each node takes into account the back propagated error and produces a revised set of weights. In this way, the network can be trained to perform its desired operation. Training may employ a supervised approach based on a set of labelled training data. Other approaches are also possible, such as a reinforcement approach wherein the network each data point is not initially labelled. The learning algorithm begins by guessing the corresponding output for each point, and is then told whether it was correct, gradually tuning the weights with each such piece of feedback. Another example is an unsupervised approach where input data points are not labelled at all and the learning algorithm is instead left to infer its own structure in the experience data.Other forms of machine learning model are also known, other than just neural networks, for example clustering algorithms, random decision forests, and support vector machines. Some machine learning models can be designed to perform causal discovery using observational data or both observational and interventional data. That is, for a set of variables (e.g. [x₁, x₂, x₃]), the model when trained can estimate a likely causal graph describing the causal relationships between these variables. E.g. in the case of three variables a simple causal graph could be x1 → x₂ → x₃, meaning that x₁ causes x₂ and x₂ causes x₃ (put another way, x₃ is an effect of x₂ and x₂ is an effect of x₁). Or as another example, x₁ → x₂ ← x₃, means that x₁ and x₃ are both causes of x2 (x2 is an effect of x1 and x3). However, in the past such models were only used for causal discovery alone, not for treatment effect estimation for decision making.Another type of model aims to perform treatment effect estimation, which commonly assumes that the causal graph is already given by the user. Such methods did not previously work with unknown causal graphs. “Deep End-to-End Causal Inference” [DECI] (Geffner et al, Microsoft Research, https://arxiv.org/pdf/2202.02195.pdf) appreciated that in many applications, users would benefit from the ability to perform treatment effect estimation for decision making with observational data only, without needing to know the causal graph. Accordingly, Geffner et al provided an integrated machine learning model that both models the causal relationships between variables and performs treatment effect estimation, by averaging over multiple possible causal graphs sampled from a distribution. This advantageously allows the method to exploit a model that has been trained for causal discovery in order to also estimate treatment effects. The method thus enables “end-to-end” causal inference. SUMMARY It is recognized herein that there is still further scope to improve on the DECI model developed by Geffner et al, or the like. Particularly, existing models do not model well time series data. The nature of cause-and-effect is such that causes in the past cause effects in the future. However, in a model based on only a static causal graph, there is nothing to stop the model learning or predicting that a cause in the future causes an effect in the past. To address such issues or similar, the present disclosure discloses a model in which, instead of a static graph that contains only edges between variables at single snapshot in time, the model samples a temporal causal graph which contains causal edges between variables at different time steps. According to one aspect disclosed herein, there is provided computer-implemented method comprising: A) selecting a selected variable from among variables of a feature vector; B) sampling a temporal causal graph from a temporal graph distribution, the temporal graph distribution specifying probabilities of directed causal edges between different ones of the variables of the feature vector at a present time step, and from one of the variables of the feature vector at a preceding time step to one of the variables of the feature vector at the present time step; and C) from among of the variables of the feature vector, identifying a present parent which is a cause of the selected variable in the present time step according to the temporal causal graph, and identifying a preceding parent which is a cause of the selected variable from the preceding time step according to the temporal causal graph. The method further comprises: D) inputting an input value of each of the identified present and preceding parent into a respective encoder, resulting in a respective embedding of each of the present and preceding parents; E) combining the embeddings of the present and preceding parents, resulting in a combined embedding; and F) inputting the combined embedding into a decoder associated with the selected variable, resulting in a reconstructed value of the selected variable. Augmenting the model to include a temporal graph will advantageously lead to more accurate predictions as it will more accurately model the causal reality of the modelled scenario. Alternatively or additionally, the model may also be used to predict the best order in which to apply a series of two or more treatments, and/or the timing with which to apply one or more treatments; and/or to predict how long a treatment will take to take effect. Another potential issue to take into account when handling time series data, is that if the modelled noise is static, then the model may not be optimally specified. For instance consider a scenario where a symptom of a modelled patient, or a sensor reading from a modelled device, remains relatively smooth during some periods, but becomes more erratic during other periods. It would be desirable to take this into account in the modelled noise. Optionally to address this, the presently disclosed model may also include a history dependent noise term, which is generated based on values of variables from past time steps. Hence in embodiments, the method may comprise generating a history dependent noise term based on embeddings of the preceding parents; and combining the history dependent noise term with the reconstructed value of the selected variable, resulting in a simulated value of the reconstructed variable. Making the simulated noise dependent on history will lead to more optimal predictions, such as estimated treatment effects, as the model will again be more likely to be a more realistic representation of the ground truth. Optionally, in some embodiments, the disclosed model may also take into account the presence of possible confounders, which have the potential to bias existing models. Consider two variables x₁ and x₂ and a model which is attempting to discover whether there is a causal edge between them, or to predict a treatment effect based on the modelled causation. A hidden confounder is a third variable which is not observed and which is a cause of both x1 and x2. This could lead to a false conclusion that x₁ is the cause of x₂ (or vice versa) when in fact there is no causal link (or a weaker causal link) and instead both x1 and x2 are an effect of a common causal variable u12 (a confounder) which is not present in the input data (i.e. not one of the variables of the input feature vector). For instance in the healthcare field x₁ could measure the presence or absence of a certain condition (e.g. disease) in a subject, whilst x2 could be a lifestyle factor such as an aspect of the subject’s diet (e.g. salt or fat intake, etc.) or whether they are a smoker, and the confounder u12 could be a measure of a socioeconomic circumstance of the patient (e.g. annual income). Ignoring the socioeconomic circumstance may give the false impression that the lifestyle factor x₂ causes the condition x1, whereas in fact the ground truth is that both x1 and x2 are effects the socioeconomic circumstance. To address such possibilities, embodiments of the present disclosure provide a machine learning model which models the possibility of hidden confounders as latent variables model. Accordingly, in embodiments B) may further comprise sampling a second causal graph from a second graph distribution, the second causal graph modelling presence of possible confounders, a confounder being an unobserved cause of both of two variables in the feature vector. In this case C) further comprises, from among of the variables of the feature vector, identifying a parent variable which is a cause of the selected variable according to the first causal graph, and which together with the selected variable forms a confounded pair having a respective confounder being a cause of both according to the second causal graph; and D) further comprises inputting the input value of the parent variable and an input value of the selected variable into an inference network, resulting in a latent value modelling the respective confounder of the confounded pair, and inputting the latent value into a second encoder, resulting in an embedding of the confounder of the confounded pair; and in E) the combining includes combining the embedding of the present and preceding parents with the embedding of the confounder of the confounded pair, thereby resulting in said combined embedding. This Summary is provided to introduce a selection of concepts in a simplified form that are further described below in the Detailed Description.This Summary is not intended to identify key features or essential features of the claimed subject matter, nor is it intended to be used to limit the scope of the claimed subject matter.Nor is the claimed subject matter limited to implementations that solve any or all of the disadvantages noted herein. BRIEF DESCRIPTION OF THE DRAWINGS To assist understanding of embodiments of the present disclosure and to show how such embodiments may be put into effect, reference is made, by way of example only, to the accompanying drawings in which: Figure 1 is a schematic block diagram of a system in accordance with embodiments disclosed herein, Figure 2 is a schematic computation diagram illustrating an example of a machine learning model, Figure 3 schematically illustrates an example of a causal graph, Figure 4 is schematic sketch of an example of a probabilistic distribution, Figure 5 schematically illustrates another example of a causal graph, Figure 6 is a schematic computation diagram illustrating a further machine learning model, Figure 7 is a schematic flowchart of a method of training a model in accordance with the present disclosure, Figure 8 is a schematic flowchart of a method of making treatment effect estimations using a trained model in accordance with embodiments disclosed herein, Figure 9 schematically illustrates a causal graph (or part thereof) including a possible confounder between a pair of variables, Figure 10 is a schematic computation diagram illustrating a machine learning model that models the possible presence of hidden confounders in accordance with embodiments disclosed herein, Figure 11 is a schematic flow chart showing the core functionality of the model of Figure 10, Figure 12 schematically illustrates an example of a temporal causal graph, Figure 13 is a schematic computation diagram illustrating a machine learning model that models temporal relationships between time series data in accordance with embodiments disclosed herein, Figure 14 is a schematic computation diagram of a history-dependent noise model in accordance with embodiments disclosed herein, Figure 15 is a schematic flow chart of a method of training the model of Figure 13, Figure 16 is a schematic flow chart of a method of using the model of Figure 13 for treatment effect estimation once trained; and Figure 17 gives additional exemplary detail expanding upon the flow chart of Figure 16. DETAILED DESCRIPTION OF EMBODIMENTS EXAMPLE SYSTEM OVERVIEW Figure 1 illustrates an example system according to embodiments of the present disclosure. The system comprises a server system 102 of a first party, a network 112, and a client computer 114 of a second party. The server system 102 and client computer 114 are both operatively coupled to the network 112 so as to be able to communicate with one another via the network 112. The network 112 may take any suitable form and may comprise one or more constituent networks, e.g. a wide area network such as the Internet or a mobile cellular network, a local wired network such as an Ethernet network, or a local wireless network such as a Wi-Fi network, etc. The server system 102 comprises processing apparatus comprising one or more processing units, and memory comprising one or more memory units. The, or each, processing unit may take any suitable form, e.g. a general purpose processor such as a CPU (central processing unit); or an accelerator processor or application specific processor such as a dedicated AI accelerator processor or a repurposed GPU (graphics processing unit), DSP (digital signal processor), or cryptoprocessor, etc. The, or each, memory unit may also take any suitable form, e.g. an EEPROM, SRAM, DRAM or solid state drive (SSD); a magnetic memory such as a magnetic disk or tape; or an optical medium such as an optical disk drive, quartz glass storage or magneto-optical memory; etc. In the case of multiple processing units and/or memory units, these may be implemented in the same physical server unit, or different server units in the same rack or different racks, or in different racks in the same data centre or different data centres at different geographical sites. In the case of multiple server units, these may be networked together using any suitable networking technology such as a server fabric, an Ethernet network, or the Internet, etc. Distributed computing techniques are, in themselves, known in the art. The memory of the server system 102 is arranged to store a machine learning (ML) model 104, a machine learning algorithm 106, training data 108, and an application programming interface (API) 110. The ML model 104, ML algorithm 106 and API 110 are arranged to run on the processing apparatus of the server system 102. The ML algorithm 106 is arranged so as, when run, to train the ML model 104 based on the training data 108. Once the model 104 is trained, the ML algorithm 106 may then estimate treatment effects based on the trained model. In some cases, after the ML model 104 been trained based on an initial portion of training data 108 and been made available for use in treatment effect estimation, training may also continue in an ongoing manner based on further training data 108 , e.g. which may be obtained after the initial training. The API 110, when run, allows the client computer 114 to submit a request for treatment effect estimation to the ML algorithm 106. The ML model 104 is a function of a plurality of variables. The request may specify a target variable to be examined, and may supply input values of one or more other variables (including intervened-on values and/or conditioned values). In response, the ML algorithm 106 may control the ML model 104 to generate samples of the target variable given the intervened-on and/or conditioned values of the one or more other variables. The API 110 returns the result of the requested causal query (the estimated treatment effect) to the client computer 114 via the network 112. In embodiments, the API may also allow the client computer to submit some or all of the training data 108 for use in the training. CAUSAL MACHINE LEARNING MODEL Figure 2 schematically illustrates an example implementation of a machine learning model 104 of a type as disclosed in “Deep End-to-End Causal Inference” (DECI) by Geffner et al. The ML model 104 comprises a respective encoder ^^ _^ ^{^} ^^{^} and a respective decoder ^^ _^ ^{^} ^^{^} for each of a plurality of variables xi, where i is an index running from 1 to D where D > 1. The set of variables xi…xD may be referred to as the input feature vector. Each variable represents a different property of a subject being modelled. The subject may for example be a real-life entity, such as a human or other living being; or an object such as a mechanical, electrical or electronic device or system, e.g. industrial machinery, a vehicle, a communication network, or a computing device etc.; or a piece of software such as a game, operating system software, communications software, networking software, or control software for controlling a vehicle or an industrial processor or machine. Each property could be an inherent property of the being or object, or an environment of the being or object which may affect the being or object or be affected by the being or object. For instance in a medical example, the variables represent different properties of a person or other living being (e.g. animal). One or more of the variables may represent a symptom experienced by the living being, e.g. whether the subject is exhibiting a certain condition such as a cough, sore throat, difficulty breathing, etc. (and perhaps a measure of degree of the condition), or a measured bodily quantity such as blood pressure, heart rate, vitamin D level, etc. One or more of the variables may represent environmental factors to which the subject is exposed, or behavioural factors of the subject, such as whether the subject lives in an area of high pollution (and perhaps a measure of the pollution level), or whether the subject is a smoker (and perhaps how many per day), etc. And/or, one or more of the variables may represent inherent properties of the subject such as a genetic factor. In the example of a device, system or software, one or more of the variables may represent an output state of the device, system or software. One or more of the variables may represent an external factor to which the device, system or software is subjected, e.g. humidity, vibration, cosmic radiation, and/or a state of one or more input signals. For instance the variables may comprise a control signal, and/or an input image or other sensor data captured from the environment. Alternatively or additionally, one or more of the variables may represent an internal state of the device, system or software, e.g. an error signal, resource usage, etc. One, more or all of the variables may be observed or observable. In some cases, one or more of the variables may be unobserved or unobservable. Each respective encoder ^^ _^ ^{^} ^^{^} is arranged to receive an input value of its respective variable xi, and to generate a respective embedding e_i (i.e. a latent representation) based on the respective input value. As will be familiar to a person skilled in the art, an embedding, or latent representation or value, is in itself a known concept. It represents in the information in the respective input variable an abstracted form, typically in a compressed form, which is learned by the respective encoder during training. The embedding may be a scalar value, or may be a vector of dimension `embedding_dim` which is greater than 1. Note that a ”value” as referred to herein could be a vector value or a scalar value. For example if one of the variables is an image (e.g. a scan of the subject), then the “value” of this vector variable is the array pixel values for the image. The different variables x_i may have a certain causal relationship between them, which may be expressed as a causal graph. A causal graph may be described as comprising a plurality of nodes and edges (note that these are not the same thing as the nodes and edges mentioned earlier in the context of a neural network). Each node represents a respective one of the variables xi in question. The edges are directional and represent causation. I.e. an edge from x_i=k to x_i=l represents that x_k causes xl (xl is an effect of xk).A simple example involving three variables is shown in Figure 3. In this example x2 causes x1, and x1 causes x3. For example x3 may represent having a respiratory virus, x₁ may represent a lung condition and x₂ may represent a genetic predisposition. Of course other possible graphs are possible, including those with more variables and other causal relationships. A given graph G may be expressed as a matrix, with two binary elements for each possible pair of variables x_i=k, x_i=l; one binary element to represent whether an edge exists between the two variables, and one to represent the direction of the edge if it exists (e.g.0 = no edge from k to l, 1 = edge from k to l, or vice versa). For example for three variables this could be written as: ^^{^ ^^ ^^1,2} ^_{^ ^^ ^^1,3} ⁾

_{^^ ^^ ^^2,3} Other equivalent representations are possible. E.g. an alternative representation of the probabilities of existence and direction of edges in a graph would be:

For any given situation, the actual causal graph may not be known. A distribution q_φ of possible graphs may be expressed in a similar format to G, but with each element comprising a parameter φ (phi, also drawn φ) representing a probability instead of a binary value. ^^{^_ ^^ ^^ ^^ ^^ ^^ ^^1,2 ^^_ ^^ ^^ ^^1,2} ^_{^ ^^ = ( ^^_ ^^ ^^ ^^ ^^ ^^ ^^1,3 ^^_ ^^ ^^ ^^1,3} ⁾ ^_{^_ ^^ ^^ ^^ ^^ ^^ ^^2,3 ^^_ ^^ ^^ ^^2,3} (Or an equivalent representation.) In other words the parameter φ_exists1,2 represents the probability that an edge between x₁ and x₂ exists; φ_dir1,2 represents the probability that the direction of the possible edge between x₁ and x₂ is directed from x₁ to x₂ (or vice versa); parameter φ_exists1,2 represents the probability that an edge between x1 and x3 exists; etc. Returning to Figure 2, the ML model 104 further comprises a selector ^^, a combiner 202 and a demultiplexer 204. It will be appreciated that these are schematic representations of functional blocks implemented in software. The selector ^^ is operable to sample a causal graph G from the distribution. This means selecting a particular graph G (with binary elements) whereby the existence and direction of the edges are determined pseudorandomly according to the corresponding probabilities in the distribution qφ. Preferably the possible graphs are constrained to being directed acyclic graphs (DAGs), for the sake of practicality and simplicity of modelling. The selector ^^ also receives a value of the index i for a selected target variable x_i. For the currently selected variable x_i (node of index i in the graph), the selector ^^ selects the respective embeddings e_Pa(i) generated by the respective encoders ^^ _^ ^{^} ^^{^} ^_{^( ^^)} of the parents Pa(i) of the node i (variable x_i) in the currently sampled graph G, and inputs these into the combiner 202. The combiner 202 combines the selected embeddings e_Pa(i) into a combined embedding e_c. In embodiments the combination is a sum. In general a sum could be a positive or negative sum (a subtraction would be a sum with negative weights). The combined (summed) representation thus has the same dimension as a single embedding, e.g. `embedding_dim`. In alternative implementations however it is not excluded that another form of combination could be used, such as a concatenation. The demultiplexer 204 also receives the index i of the currently selected variable x_i, and supplies the combined embedding e_c into the input of the decoder ^^ _^ ^{^} ^^{^} associated with the currently selected variable x_i. This generates a value of a respective noiseless reconstructed version x_i’ of the respective variable xi based on the combined embedding ec. The encoders ^^ _^ ^{^} ^^{^} =_{1… ^^} and decoders ^^ _^ ^{^} ^^{^} =_{1… ^^} are constituent machine learning models comprised by the overall ML model 104. They are each parameterized by respective sets of parameters θ which are tuned during learning. In embodiments, each of the encoders ^^ _^ ^{^} ^^{^} =_{1… ^^} and decoders ^^ _^ ^{^} ^^{^} =_{1… ^^} is a respective neural network (in which case their parameters may be referred to as weights). However it is not excluded that some or all of them could instead be implemented with another form of constituent machine learning model. TRAINING Figure 7 schematically represents a method of training the DECI type model 104 based on a training data 108. The training data 108 comprises a plurality of training data points. Each data point [x1 … xD] comprises a set of input values, one respective value for each of the variables xi. The method processes each of the training data points, at least in an initial portion of the training data 108. At step S10 the method takes a new data point from memory, and at step S20 the index i is set to a first value (e.g. i=1) to specify a first target variable xi from among the set of variables in the data point. At step S30, this causes the selector ^^ to select the target variable xi with the currently set value of the index i as the variable to be processed. At step S40 the selector ^^ samples a random graph G from the distribution q_φ. At step S50, the selector ^^ selects the parents Pa(i) of the target variable xi (node i in the graph) and supplies the respective embeddings ePa(i) from the encoders ^^ _^ ^{^} ^^{^} ^_{^( ^^)} of the selected parents into the combiner 202. At step S60, the combiner 202 combines (e.g. sums) the embeddings of the selected parents Pa(i) into the combined embedding e_c, and the demultiplexer 204 selects to supply the combined embedding e_c into the decoder ^^ _^ ^{^} ^^{^} of the target variable xi. The respective decoder ^^ _^ ^{^} ^^{^} is thus caused to generate a noiseless reconstruction xi’ of the selected target variable xi. This could equally be expressed as:

where ^∑ _{^^∈ ^^ ^^( ^^)} is another way of writing the operations of the selector ^^ and combiner 202.Preferably, G is a DAG (directed and acyclic) so that it is valid to generate the value of any node in this way. The difference x_i-x_i’ may be referred to as the residual. The process is repeated for each variable xi in the set of variables, i.e. for i=1…D, thus creating a full reconstructed set of values [x₁’ … x_D’] for the data point in question. This is represented schematically by steps S70 and S80 in Figure 7, looping back to S30. However note that this the “loop” may be a somewhat schematic representation. In practice some or all of the different variables xi for a given data point could be processed in parallel. At step S90 the ML algorithm 106 applies a training function which updates the parameters (e.g. weights) θ of the encoders ^^ _^ ^{^} ^^{^} =_{1… ^^} and decoders ^^ _^ ^{^} ^^{^} =_{1… ^^} . Simultaneously it also updates the parameters φ of the distribution q_φ of possible graphs. The training function attempts to update the parameters θ, φ in such a way as to reduce a measure of overall difference between the set of input values of the set of input variables [x1 … xD] and the reconstructed version of the set variables [x₁’ … x_D’]. For example in embodiments the training function is an evidence lower bound (ELBO) function. Training techniques of this kind are, in themselves, known in the art, e.g. based on stochastic back propagation and gradient descent. During the training phase, x_i is an observed data value that is attempted to be reconstructed. The residual x_i—x_i’ between the reconstruction and the observation may be used to compute the ELBO objective. In embodiments, the model 104 may be a probabilistic model that it assigns a joint probability (likelihood) p (X₁,…,X_D) to all the variables X₁,…,X_D; where X_i represents the noiseless reconstruction x_i combined with a random noise term z_i, e.g. additively (X_i = x_i’ + z_i), and zi may be sampled from a random distribution, zi ~ p(zi). In such embodiments, then in order to train the model 104, one first collects a real-world data point with D-dimensional features (x1, x2,…,xD), and then maximizes the likelihood of it being generated by the model. That is to say, the training operates to try to maximize p(X1 = x1, X2 = x2,…,XD = xD). This is done by adjusting the parameters (e.g. weights) θ of the model 104. The ELBO mentioned previously is an example of a numerical approximation to the quantity p(X₁ = x₁, X₂ = x₂,…,X_D = x_D). ELBO is typically much easier to compute than p(X1 = x1, X2 = x2,…,XD = xD). The above-described process is repeated for each of the data points in the training data 108 (or at least an initial portion of the training data). This is represented schematically in Figure 7 by step S100 and the loop back to S10. Though again, this form of illustration may be somewhat schematized and in some implementations, batches of data points could be reconstructed in parallel, and the model parameters updated based on the batches. At the beginning of the overall training method the graph distribution q_φ may start with some predetermined set of probabilities, e.g. all elements 0.5 for the existence and direction of each possible edge, or using some prior domain knowledge to inform which edges are possible or impossible, or more or less likely. Then, as the model is trained with more and more data points, the probabilities (i.e. the parameters φ of qφ) are gradually learned in parallel with the parameters θ (e.g. weights) of the encoders and decoders ^^ _^ ^{^} ^^{^} =_{1,…, ^^} , ^^ _^ ^{^} ^^{^} =_{1,… ^^} . TREATEMENT EFFECT ESTIMATION Once the model 104 has been trained sufficiently based on all the data points in the training data 108 (or at least an initial portion of the initial training data), then at step S110 the trained ML model 104 is made available to be used for treatment effect estimation (i.e., answering causal queries). In embodiments, this may comprise making the model 104 available to via the API 110 to estimate treatment effects / answer causal queries requested by the client computer 114. A method of using the ML model 104 to estimate a treatment effect is shown schematically in Figure 8. At step T10 the index i is set to that of a target variable xi whose treatment effect is to be estimated. In addition, the input values of one or more ”intervened-on” variables are also set to their known values. The intervened-on variables are variables whose values are set to some specified value, to represent that the property that they represent has been controlled (the treatment, i.e. an intervention on the modelled property). An ”intervened-on” variable could also be referred to as a treated variable or controlled variable. For example, in the medical application, the intervened- on variable(s) may represent one or more interventions performed on the subject, and the target variable may represent a possible symptom or condition of the subject (e.g. the presence of a certain disease). Or in the case where the subject is a device or software, the intervened-on variable(s) may represent one or more states that are set to defined values, and the target variable may represent a condition or state of the device or software that is being diagnosed. At step T20, the selector ^^ samples a graph G pseudorandomly from the distribution q_φ according to the probabilities φ. In other words the chance of there existing an edge in given direction between each pair of nodes i=k, i=l in the sampled graph G is determined by the pair of parameters φ _exists(k,l); φ_dir(k, l). Although optionally, some of the learned probabilities could be overridden to predetermined values based on prior knowledge (e.g. in some scenarios a given causality could be ruled out – probability set to zero – or may be known to be unlikely, either based on a priori or empirical knowledge). At step T30 the selector ^^ selects the parents Pa(i) of the target variable xi in the sampled graph G, and supplies the respective embeddings ePa(i) from the encoders ^^ _^ ^{^^} ^_{^^( ^^)} of the selected parents into the combiner 202 to be combined (e.g. summed) into the combined embedding ec. The demultiplexer 204 selects to pass the combined embedding ec into the decoder ^^ _^ ^{^} ^^{^} of the target variable xi, thus causing it to generate a noiseless reconstruction of xi’. Thus the selection ^^ selects the parents of node i in the graph G, so it depends on both i and G. The index is set to a particular target node i that one is trying to estimate (as well as selecting a graph G from q_φ), and the model 104 outputs a noiseless reconstruction of x_i, denoted by x_i’. During training, such reconstruction was formed of each node given the actual data for its parent nodes. Now at test time, the model 104 is used in a slightly different way for treatment effect estimation.Here, the target variable x_i which is now treated as unknown, and the goal is to generate simulated values of the target variable, given other interventional variables. In some embodiments the simulated value could just be taken as the noiseless reconstructed value x_i. However for some purposes, the noiseless reconstruction x_i’ of the target variable x_i may not be considered preferable, and instead the full interventional distribution may be taken into account. This can characterized by the following sources of uncertainties: 1), different realizations of the graph G that could happen to have been selected during sampling (even an unlikely one), which can be modelled by simulating different graphs from q_φ.; and 2), the randomness of the residual noise random variables zi. The latter can be taken into account by the following simulation equation:

In other words, during the treatment effect estimation phase, X_i is a simulated value that is obtained through random sampling of zi, wherein Xi = xi’ + zi (except for xi=T, the intervened-on variable). More generally in other possible implementations, the noise z_i could be combined with the noiseless reconstruction x_i’ in other ways in order to obtain Xi, not necessarily additively. In the following description, reference may be made to Xi as the simulated value of target variable xi. However generally the simulated value could simply be the noiseless reconstructed value xi’, or X_i (which takes into account noise, e.g. X_i = x_i’ + z_i), or a value which is simulated based on the reconstructed value xi’ in any other way. By whatever means the simulated value is determined, according to the present disclosure, an average is taken over multiple sampled graphs.In embodiments, an average is taken over multiple sampled graphs and residual noise variables zi. In other words, for a given target variable xi to be estimated, multiple values of variable Xi are simulated based on different respective sampled graphs G, each time sampled randomly from both q_φ and z_i This is represented schematically in Figure 8 by the loop from step T50 back to T20, to represent that multiple simulated values of xi are determined (for a given target variable xi), each determined based on a different randomly sampled graph G (each time sampled from the distribution qφ) and optionally noises ^^ _^^ . Again the illustration as a loop may be somewhat schematized, and in embodiments some or all of the different simulated values for a given target variable x_i, based on the different sampled graphs and sampled residual noises z, may in fact be computed in parallel. In embodiments the average could be taken as a simple mean, median or mode of the different values of the simulated values x_i’ or X_i’ of the variable x_i (as simulated with the different sampled graphs and optionally residual noises). In embodiments, such averaging is based on estimating an expectation of the probabilistic distribution of the simulated values. Figure 4 schematically illustrates the idea of a probabilistic distribution p(xi). The horizontal axis represents the value of the target variable being simulated, x_i, and the vertical axis represents the probability that the variable takes that value. If the distribution is modelled as having a predetermined form, it may be described by one or more parameters, e.g. a mean μ and standard deviation σ or variance σ² in the case of a Gaussian. Other more complex distributions are also possible, which may be parameterized by more than two parameters, e.g. a spline function. As a convenient shorthand the probabilistic distribution may be referred to as a function the target variable of x_i, i.e. p(xi), in the sense that it models the distribution of the target variable x_i, but in fact it will be appreciated that the distribution is in fact determined based on the corresponding simulated values xi’ or Xi. Based on the trained model 104, it is possible using statistical methods to estimate an expectation E[p(x_Y)|do(x_T=val1)] of the probabilistic distribution p of a particular target variable x_i=Y given the intervened value val1 of another, treatment/intervention variable xi=T. In embodiments, the target variable xY may model an outcome of a treatment modelled by xT. Note that “treatment” as used most broadly herein does not necessarily limit to a medical treatment or a treatment of a living being, though those are certainly possible use cases. In other examples, the treatment may comprise applying a signal, repair, debugging action or upgrade to an electronic, electrical or mechanical device or system, or software, where the effect may be some state of the device, system or software which is to be improved by the system. In embodiments, the actual real-world treatment may be applied in dependence on the estimation (e.g. expectation) of the effect of the modelled treatment, for example on condition that the treatment estimation (e.g. expectation) is above or below a specified threshold or within a specified range. One suitable measure of expectation may be referred to as the average treatment effect (ATE). This may be expressed as: E[p(x_Y)|do(x_T=val1)] - E[p(x_i=Y)|(do(x_T=val2)] or E[p(xY)|do(xT=val1)] - E[p(xi=Y)] which could also be expressed E[p(xY)|do(xT=val1)] - E[p(xi=Y)|do(xT=mean(xT))], where val1 is some treatment, val2 is some other treatment, and “do” represents applying the treatment. In other words, the ATE is the difference between: a) the expectation of the distribution of the effect xY given the value val1 of a treatment xT and b) the expectation of the distribution of the effect x_Y given the value val2 of a treatment x_T, or the difference between a) the expectation of the distribution of the effect xY given the value val1 of a treatment xT and b) the expectation of the distribution of the effect xY without applying a value of the treatment xT. Note that in variants the equality x_T=val1 or x_T=val2 could be replaced with another form of expression such as a “greater than”, “less than” or range type expression. An estimation of E[p(xY)|do(xT=treatment) may be determined based on the modelled distributions of noise, together with the learned graph and arrow functions. This may be done using the noise distribution at node x_Y, and using noise samples from other upstream variables, plus knowledge of the graph and its behaviour under the intervention do(xT=treatment), in order to thus make the ATE estimation calculation. The expectation E[p(xY)|do(xT=val)] of a target variable x_Y given a treatment x_T=val may be computed as:

where z_Y is a noise term. E is an average over everything that is random in the equation. Consider for example a simple graph: x1→ x2 ←x3 where x1 is the treatment xT, and x2 is the target variable xY. In the computation of E, x1 is set to its known value. E may then be written: ^^[ ^^ _^ ^{^} ^^{^}( ^^₁ ^{^^}( ^^₁ = ^^ ^^ ^^) + ^^₃ ^{^^}( ^^₃)) + ^^ _^^] To compute E, one example method is simply to take the mean of the different sampled values xY’ or XY of the target variable xY, based on the different sampled graphs, optionally also including some random noise (e.g. additive noise) in each sample. Another option is to fit the sampled values of x_Y to a predetermined form of distribution, such shown as in Figure 4, e.g. a Gaussian, normal or spline function (again optionally also including random noise). The average may then be determined as the average of the fitted distribution. Note that during the determination of the multiple instances of x_i’ it is possible, especially for some distributions, that the same graph G ends up being sampled multiple times, or even every time.However the graph is still being freshly sampled each time and even if the sampled graphs turn out to be the same, they may still be described as different instances of the sampled graph, in the sense that they are individually sampled anew each time. In the case where the multiple graphs happen to be the same, then the averaging only averages the noise z. In some cases, if the treated (i.e. controlled) variable has one or more parents, the graph G may be mutilated: G->G_do(xT=val). E.g. consider the graph: x4→ x1→ x2 ←x3 where again x₁ is the treatment x_T, and x₂ is the target variable x_Y. In the computation of E, x₁ is set to its known value. Therefore x₄ has no effect on the outcome of x₂. So in the determination of the expectation, the graph is mutilated to remove the node x4 and the edge from x4→ x1. In other words, fixing the value of the known (controlled) variable x1 means that any effect of the edge from the parent of the known variable x₁. Another type of average which may be determined according to embodiments disclosed herein may be referred to as the conditional ATE (CATE). This is the difference between: a) the expectation of the target variable x_Y given the treatment x_T=val1 conditional on an observation γ of at least one other of the variables x_C in the set, and b) the expectation of the target variable x_Y without the treatment (either with a different treatment val2 or no treatment) but still conditional on the same observation of x_C. That is, E[p(x_Y)|do(x_T=val1), x_C=γ] - E[p(x_Y)|do(x_T=val2), x_C=γ] or: E[p(xY)|do(xT=val1), xC=γ] - E[p(xY)|,xC=γ] which could also be expressed E[p(x_Y)|do(x_T=val1), x_C=γ] - E[p(x_i=Y)|do(x_T=mean(x_T)) x_C=γ]. Note that in variants the equality x_T=val1, x_T=val2, or x_C=γ could be replaced with another form of expression such as a “greater than”, “less than” or range type expression. Figure 5 illustrates by way of example why estimating the conditional treatment effect is not necessarily straightforward. In this example x₃ is the target x_Y whose treatment effect will be estimated and x4 is the treatment xT, where the treatment is the cause of the target effect. In addition, there is another, unobserved cause x₁ of the target effect, and another observable effect x₂ of the unobserved cause x₁. The variable x₂ is to be the observed condition x_C. For instance, in a medical example the target effect x3 could be some condition or symptom of the subject (e.g. a respiratory problem), the treatment x4 could be a possible medical intervention (e.g. taking some drug or vitamin), the unobserved cause x₁ may be a genetic factor, and the other observable cause x2 may be some observable physical quality of the subject’s body (e.g. body mass index). In general the unobserved cause could be unobservable, or merely unobserved. The observable condition x₂ contains information about the unobserved cause x₁, which in turn can reveal information about the desired effect x3 (=xY). For example if it is known that an athlete is in the Olympics and that they are a footballer, this reduces the probability that they are also a rower. In fact for any two variables that are effects of a common cause and whose intersection is not 100%, knowing something about one gives information about the other. E.g. if it is known that a subject has lung cancer and that they were exposed to a carcinogenic chemical other than tobacco, this reduces the probability that they were are smoker. However the causal direction is from x₁ → x₂, and the model 104 of Figure 2 is only configured to learn effects of causes in the direction from cause to effect – it is not configured to learn inferences of effect from cause. I.e. it is not configured to “go against the arrows” in the figure (the directional causal edges). To address this, as illustrated in Figure 6, in embodiments the ML model 104 may be adapted to include at least one inference network h disposed between at least one observable condition xc (x2 in the example) and at least one unobservable potential cause (x₁ in the example) of the condition x_c. As illustrated in Figure 6, in some such embodiments, the inference network h (or individual such networks) may be disposed between the unobserved cause and multiple potential effects (up to all the other variables). This will allow the model to learn which variable(s) may be an effect of the unobserved cause, if relationship is not prior knowledge. The inference network(s) h may be trained at the training stage simultaneously along with the encoders g^e and decoders g^d and the parameters of the graph distribution qφ, or alternatively after the rest of the model (see below). In embodiments the inference network h may comprise a neural network, in which case training the inference network comprises tuning the weights of the inference network. Alternatively the use of other forms of machine learning is not excluded for the inference network. The inclusion of the inference model makes it possible to estimate a conditional expectation such as E[xY|do(xT=val1), xC=γ]. To estimate the conditional expectation, then initially during the estimation of the target variable xY, the conditional variable xc is not fixed. In other words, the method proceeds as described above with respect to ATE, to obtain multiple different samples of x_Y based on multiple respective sampled graphs. At the same time, respective simulated samples x_c of the conditional variable are also obtained in the same way based on the respective sampled graphs. This gives a two dimensional set of samples (xY, xc), which could be described as a plot or “scatter plot”.Then a predetermined form of function is fitted to the 2D set of samples (x_Y, x_c), such as a straight line, a curve, or a probabilistic distribution. After that, xc is set to its observed value in the fitted function, and a corresponding value of xY is read out from the fitted function. This is taken as the conditional expectation of x_Y given x_c. At least two alternative variants to computing CATE may be employed, depending on implementation. Variant I: estimate CATE using the same approach as used to estimate ATE, but performing a re- weighing of the terms inside of the expectations such that the condition (that gives CATE its name) is satisfied. This type of approach is known as an importance sampling technique. The weights of the different samples used to compute the expectation are provided by an inference network, which is trained together with the rest of the ML model 104. Variant II: after the model 104 of Figure 2 has been trained and a specific CATE query is received (e.g. via the API 110), the inference network h is trained to estimate the effect variable from the conditioning variable. To train this model, data simulated from the trained model of Figure 2 is used, while applying some treatment specified in the query. Then the conditional average treatment effect is estimated by inputting the relevant value of the conditioning variable into the inference network h. It returns a distribution over effects from which the expected effect can be computed. The main difference between the two approaches is that the first solely uses components learnt from the observed data during model training, while the second requires learning a new network after the model 104 has been trained. The reasoning for proposing both methods is that there are specifics settings where one or the other are more computationally efficient. Note that the disclosed methods are not limited to the controlled (i.e. treated) variable xT being a direct parent of the target variable x_Y (i.e. the variable whose treatment effect being estimated). The treated variable or the effect variable can be any arbitrary variable in the set xi=1…D, e.g. a grandparent of the target variable xY. In embodiments, the simulation of the target variable takes into account a potential effect of all causal variables across the sampled graph. An example implementation of this is as follows. This may be used in conjunction with any of the ways of averaging discussed above, or others. The method of estimating the target variable xY (e.g. the treatment effect) may comprise an inner and an outer loop. In the inner loop, the method loops through i=1…D, so as to generate the simulated value X_i (or just x_i’) of each of the input variables x_i (even those that are not the target variable), except skipping over the one or more controlled (treated or known) variables x_T which are set to their known, fixed values (any edges into those nodes have been “mutilated” away, i.e. removed from the sampled graph). Then, proceeding to the next round (iteration) of the outer loop, the simulated values X_i of the non-controlled variables are fed back to the respective inputs of the model 104. In other words (other than for the one or more controlled variables xT), the simulated values Xi from the previous round or cycle (iteration) of the outer loop become the input values xi of the current iteration of the outer loop to generate an updated set of values for the simulated variables xi. This may be repeated one or more further times, and the simulated values will start to converge (i.e. the difference between the input layer and output layer of the model 104 will get smaller each time). If noise is included the noise is frozen throughout a given inner loop, then re- sampled each outer loop. The total number of iterations of the outer loop may be predetermined, or the outer loop may be iterated until some convergence criterion is met. In embodiments the outer loop is iterated at least D-1 times, which guarantees convergence without needing to evaluate a convergence criterion. This method advantageously allows causal effects to propagate throughout the graph. For example if x1 causes x2 and x2 causes x3, and an intervention is performed on x1 then the outer loop will be run at least two times to propagate the effect through to x3. However the method of using an inner and outer loop is not essential. An alternative would be to perform a topological sort of the nodes and propagate effects through in a hierarchical fashion starting from the greatest grandparents or “source nodes” (those which are only causes and not effects). Or as another alternative, though less preferred, it is not necessary to propagate effects through multiple generations and instead only the effects of immediate parents of the target variable may be taken into account. OPTIONAL EXTENSION TO INCLUDE LATENT CONFOUNDERS The following describes an optional extension to the model to accommodate the possible presence of “confounders”. Figure 9 schematically illustrates the issue of the presence of possible “hidden confounders”, which may also be referred to just as confounders. A confounder is a variable that is both A) unobserved (i.e. hidden, so not part of the input feature vector); and B) a cause of at least two variables that are among the variables of the input feature vector [x₁….x_D]. These two variables may be referred to herein as the “confounded” variables or a “confounded pair” of the respective confounder. In the example of Figure 9 the confounded variables are labelled x1 and x2 and the confounder between them is labelled u12. The confounder could be unobservable or merely unobserved. Of course the causal graph may be more complex than shown in Figure 9, and confounders may exist between more than one pair of variables of the input feature vector. An issue recognized herein is that existing machine learning models do not take into account the presence of possible confounders, which may lead to an erroneous bias in the training of such models or predictions made by such models. Consider for example a scenario where x1 is a variable measuring the presence or absence of a certain condition (e.g. disease) in a subject (e.g. patient), and x2 is a variable measuring a property of the subject which may or may not affect the possible presence of the condition. For instance x2 could measure a lifestyle factor such as an aspect of the subject’s diet (e.g. salt or fat intake, etc.), whether they are a smoker, a dosage of a certain medication they are taking, or an environmental factor such as living in an area of high pollution or near an electricity pylon. Say then that the model is being used to try to learn whether there is a causal link between x₁ and x₂, or to predict what would be the effect on x₁ of treating x₂. However there may be an unobserved common cause of the two variables, i.e. a confounder. For instance, the confounder u12 could represent a factor in the socioeconomic circumstances of the patient (e.g. annual income, education or family circumstances). Ignoring the confounder (e.g. socioeconomic circumstance) may give the false impression that the lifestyle factor x₂ causes the condition x1, whereas in fact the ground truth is that both x1 and x2 are effects of the socioeconomic circumstance, in which case x2 may not actually be a cause of x1 or may only be a weaker cause than it would otherwise appear. The presently disclosed extension – which in embodiments may be applied as an extension to the previous work by Geffner et al summarised in the Background section – introduces the modelling of possible hidden confounders into a causal machine learning model by introducing a second causal graph sampled from a second graph distribution. The extended model works in a similar manner to the model 104 described in relation to Figure 2 and like elements function in the same manner except where stated otherwise. Figure 10 schematically illustrates the extension to the model on the encoder side. The decoder side works the same as described in relation to Figure 2. As illustrated in Figure 10, as in the base model 104 of Figure 2, the extended model 104’ comprises a respective first encoder g^e _i for each variable x_i (i.e. each feature) of the input feature vector [x1…xD], each arranged to encode the value of its respective variable xi into a respective embedding ei. The extended model 104’ additionally comprises an inference network H, and a respective second encoder g^e _ij for each of a plurality of pairs of input variables x_i, x_j in the input feature vector – preferably one for each possible pair (i.e. every combination of two variables from the input vector). For each of the plurality of pairs of variables xi, xj, the inference network H is arranged to encode the respective values of one or more of the input variables x1…xD (e.g. in training, at least the two respective values of the pair together) into a respective latent value u_ij for the respective pair. This latent value represents a modelled value of the (possible) confounder that potentially exists between the respective pair of input variables xi, xj. The inference network H may also be a function of a pseudorandom noise term εi. These noises will be added as part of the input. That is, u = H(x, ε), as shown in Figure 10. In embodiments, the inference network H is implemented as one common inference network that encodes all the variables x1…xD together into the latent confounder value uij for a given pair. Alternatively it is not excluded that the inference network H could be structured as a respective constituent inference network Hij for each of the plurality of pairs of input variables input variables xi, xj. Other internal structures of H are also possible. The extended model 104’ further comprises a respective second encoder g^e _ij for each respective latent value u_ij (corresponding to each respective pair of input variables x_ij). Each respective second encoder g^eij is arranged to encode the respective latent value uij of the respective confounder into a respective embedding eij. The extended model 104’ comprises a selector ^^’ analogous to the selector ^^ described in relation to the base model of Figure 2, but with additional functionality. As in Figure 2, a causal graph G_φ is sampled from the graph distribution q_φ, and a certain selected variable x_i is selected; and the selector ^^’ then selects the parent variables Pai that are parents of the selected variable xi according to sampled graph Gφ, and inputs these into the combiner 202 (e.g. adder or concatenator). The graph distribution used for this is the same as the graph distribution q_φ described previously with respect to Figure 2. I.e. it comprises matrix whose elements represent the probabilities of directed causal edges existing between the different possible pairs of variables in the input vector [x1 … xD]. In the context of the extended model 104’, this may also be described as the first graph distribution, or the directed graph distribution; and the graph sampled therefrom may be referred to as the first causal graph or the directed graph. In other words: qφ = q1 = qdirected G_φ = G₁ = G_directed However, in addition, the extended model 104’ also samples a second causal graph G₂ from a second graph distribution q2. These may also be referred to herein as the bidirected graph distribution and bidirected graph: q₂ = q_bidirected G2 = Gbidirected The second graph distribution q2 is a graph distribution representing the probabilities that confounders exist between pairs of variables. It may be expressed as a matrix in which the different elements represent the probabilities that there exists a confounder between different pairs of variables x_i, x_j, preferably having an entry for each possible pair of variables in the input feature vector [x_i…x_D] (i.e. every combination of two variables from the input vector). So for a simple three variable input vector (D=3), the second graph distribution would have three entries: Ψ_exists12 Ψ_exists13 Ψ_exists23 where Ψ_exists(i,j) represents the probability that a confounder exists between variables with indices i and j. More generally the second graph distribution q₂ may be expressed as:

Other equivalent representations may also be possible. To sample an actual graph G₂ from the second graph distribution q₂, the presence or absence of each edge in the graph G2 is determined pseudorandomly according to the probability specified in the corresponding element of the distribution q2. E.g. if ψ_exists1,1 = 0.7, then the element (1,2) in G2 – representing whether or not a confounder u₁₂ exists between variables x₁ and x₂ – has a 70% chance of being 1 (meaning confounder present) and a 30% chance of being 0 (meaning no confounder present). Again, other equivalent mathematical expressions of the same functionality are also possible. The value of the confounder, if it exists in the sampled graph G₂, taken from the output of the inference network H. In embodiments H may be configured to model a probabilistic distribution, which may also be described as a third distribution q3, from which the values of uij are sampled. Based on the sampled graphs G₁ and G₂, in addition to the embeddings e_i of the parents Pa(i) of the selected variable xi determined from the first causal graph G1 (=Gφ), as sampled from the first graph distribution q1 (=qφ), the selector ^^’ also selects the embeddings eij of any confounded pairs that include the selected variable xi - i.e. of any confounders uij that affect the selected variable xi - as determined from the second causal graph G₂ sampled from the second graph distribution q₂. These embeddings eij of the selected confounders uij are input into the combiner 202 (e.g. adder or concatenator) along with the embeddings ei of the parents Pa(i) of the selected variable xi. The combiner 202 then combines (i.e. sums or concatenates, depending on implementation) all of the embeddings e_i of the selected parents and the embeddings e_ij of the selected confounders together into a combined embedding eC. As shown in Figure 2, the combined embedding eC is then input into the decoder g^di associated with the selected variable x_i, which decodes the combined embedding to produce a respective reconstructed value xi’ of the input variable xi. Optionally the pseudorandom noise zi may be added (or included in some other, non-additive way) to create the noisy reconstructed value X_i. Either may be used as the simulated value of x_i. This side of the model 104’ (the decoder side) works the same as in the base model 104, as described previously with respect to Figure 2. It will be appreciated that the inference network H, second encoders g^eij, and selector ^^’ may be implemented as modules of software in a similar manner to the other components as described in relation to the model of Figure 2. As with the first encoders g^e _i and decoders g^d _i, each of the inference network H and second encoders g^e _ij may be implemented as a respective neural network, or alternatively other types of constituent ML models such as random forests or clustering algorithms, or any other form of parametric or non-parametric function, are not excluded. The core functionality of the extended model 104’ is represented in the flow chart of Figure 11. By core functionality here, it is meant the mechanism of the model for sampling a causal graph that includes latent confounders, identifying parents of a selected variable according to the sampled graph, including any confounders, and then combining and decoding the embeddings of the identified parents, including the embeddings of any confounders. This mechanism may be used in both training and in subsequent treatment effect estimation. At step U10 the method selects one of the variables x_i from among the features of the input feature vector [x₁…x_d]. When used in training this variable x_i is one of the variables to be reconstructed in order to determine a loss effect. When used in treatment effect estimation the selected variable x_i will be the target variable, and the intervened-on variable (i.e. variable to be treated) will also be selected separately (though may be null). At step U20 the method samples a first causal graph G1 from the first graph distribution q1 and samples a second causal graph G2 from the second graph distribution q2. At step U30 the method determines which of the other variables [x1 … xi-1, x_i+1 … x_D] in the feature vector are parents Pa(i) of the selected variable x_i in the sampled first graph G1. Note this may involve mutilating the graph as described previously. The method also determines which of these parents Pa(i) share a confounder uij with the selected variable xi in the sampled second graph G₂. In general this could be none, one, some or all of the parents Pa(i) identified from the first graph G1. At step U40 the method generates the reconstructed value xi’ of the selected variable xi. This is done by inputting the input values of each of the selected parents Pa(i) (as determined from the first graph G₁) into the respective first encoders g^e _i in order to generate their respective embeddings ei; and inputting at least the input value of the parent in each selected confounded pair xi, xj (as determined from the second graph G2) into the inference network H to generate the respective latent value u_ij, and inputting this into the respective second encoder g^eij for the respective confounder in order to generate the respective embedding eij for each respective pair of confounded variables. In some circumstances, the values input into the inference network H to generate the latent value u_ij for a confounded pair i,j (and in turn the respective embedding e_ij) may also comprise one or more additional input values in addition to the parent of the respective pair. E.g. in training these one or more additional input values preferably comprise the observed input value of the selected variable as well. And in training or inference (e.g. treatment effect estimation), the one or more additional input values may comprise input values of one or more other of the variables of the feature vector {x1…xD} (other than the selected variable xi or the respective parent xj of the pair i,j), as one or more of these variables may potentially also comprise information about the possible confounder. In general, H can generate u_ij for the selected variable xi based on what is known about any or all of the other variables. The more is known, the better that H will predict uij. One (of many) ways to it is to take the average: H(xj) = sum over xi {H(xi, x_j) * p(x_j)}. So H can generate u_ij from only a specified value of the parent x_j. H can even generate u_ij given nothing as input, meaning H is just generating u_ij randomly from a distribution learned at training. When estimating average treatment effects, this is acceptable, since in the one will marginalize out all the variables x in the entire population anyway. Continuing step U40, all the generated embeddings e_i; e_i,j for the selected parents and confounders are input into the combiner 202 to be combined (e.g. summed or concatenated) into the combined embedding eC, which is then decoded by the respective decoder g^di for the selected variable xi to produce the respective reconstructed value x_i’. Optionally a noise term z_i may be added to produce the simulated value Xi. Alternatively the noiseless reconstructed value xi’ could be used as the simulated value. It will be appreciated that Figure 11 is somewhat schematized and all the steps of Figure 11 do not necessarily have to be performed linearly. E.g. the method could identify the parents Pa(i) from the first graph and begin generating the parents’ embeddings ei before it has begun identifying all of the confounded pairs xi, xj. Or the method could begin generating embeddings for some parents while still selecting other parents from the graph G₁, etc. Training of the model 104’ may be performed in an analogous manner to that already described in relation to Figures 2 and 7, by replacing steps S30-S70 in Figure 7 with steps U10-U40 from Figure 11. In other words, training begins with a given training data point comprising a set of respective values of at least some the variables (features) of the feature vector xi = [x1 … xD]. Preferably this will be values of all the variables of the feature vector if available, but the function H can take any number of observations as input, so if one or more values are unobserved then they are simply not passed through H. The method then cycles through the indices of each variable i = 1…D observed in the feature vector, and for each performs the method of Figure 11 with the current variable as the selected variable in order to produce a respective simulated value xi’ or Xi. Over the cycle, the method will thus produce a simulated value xi’ or Xi for each variable in the feature vector (or at least those that are observed), which together thus result in a simulated feature vector x_i’ = [x₁’ … x_D’] or X_i = [X₁ … X_D]. The method then determines a measure of the discrepancy between the input feature vector xi and the simulated feature vector xi’ or Xi for the given training data point, and based on this tunes the parameters (e.g. weights) of all the constituent models (e.g. neural networks) of the ML model 104’, including the first encoders g^e _i, second encoders g^eij, inference network H, and decoders g^di, as well as the first and second graph distributions q1 and q2. E.g. this may be done using an ELBO function, or any other suitable learning function. Algorithms for comparing an input feature vector with a simulated version in order to tune the parameters of a ML model are, in themselves, known in the art. The method may be repeated over multiple training data points in order to refine the training. Once the model 104’ has been trained to at least some extent, it may be used to estimate whether a likely confounder exists between a given target pair of variables x_i, x_j (i.e. a pair of variables of interest). This may be done by accessing the trained second graph distribution q2 from memory and reading out the value of the probability ψ_exists(i, j). If the probability is high (e.g. above a threshold) it may be determined that an unobserved confounder is likely to exist between the two variables xi and xj, whereas conversely, if the probability is low (e.g. below a threshold) it may be determined that such a confounder is not likely to exists between the variables in question. In dependence on this, a decision may be made as to whether or not to action a treatment on a first, intervened-on one of the two variables xj as a means of affecting a second, targeted one of the two variables xi. Consider a scenario where xj is a apparently cause of xi according to the first graph distribution q₁, because the relevant elements – e.g. φ_exists(I,k), φ_dir(i,j) – of the first graph distribution q₁ have been read out from memory and indicate a high probability (e.g. above a threshold probability) of a directed edge existing from xj to xi. Thus from the first graph distribution q1 alone, it may appear that treating xj is an effective way of affecting xi. However, if ψ_exists(i, j) in the second graph distribution q₂ also indicates a high probability (e.g. above a threshold probability) of a hidden confounder existing between xj and xi, then it may be determined that in fact xj should not be treated as a means of affecting target variable xi. The trained model 104’ may also be used to estimate the value of such a confounder. This may be done by reading out the value of uij from the output of the inference network after a suitable amount of training. In embodiments the inference network H may model a probabilistic distribution, which could also be described as q₃, from which the value u_ij is sampled. The decision as to whether to action the treatment may also be dependent on the value of the estimated confounder. E.g. if a confounder is determined likely to be present but having a weak effect, it may be decided still to action the treatment of xj as a means of affecting target variable xi. In alternative or additional applications, the trained model 104’ may be used to perform treatment effect estimation in an analogous manner to that already described in relation to Figures 2 and 8, by replacing steps T10-T40 in Figure 8 with steps U10-U40 from Figure 11. In other words, the target variable is selected as the selected variable xi, and input values of one or more other, intervened-on (i.e. controlled) variables xj is/are set to its/their intended value (i.e. the proposed value for treatment). As discussed earlier the inference network H can take any number of observations, so if any variables are not intervened-on or observed, their values are simply not passed through H. The method of Figure 11 is then used to generate a simulated value xi’ or Xi of the target variable x_i. In embodiments this may be repeated over multiple different sampled instances of the first and second causal graphs G1, G2; i.e. repeating the steps of Figure 11 but newly re-sampling the first causal graph G1 from the first graph distribution q1 and the re-sampling the second causal graph G₂ from the second graph distribution q₂ afresh each time (while keeping the same selection of the target variable and the same fixed values of the intervened on variable or variables xj). The estimated effect of the treatment may then be determined by averaging over these multiple rounds of graph sampling. E.g. this may comprise determining the average treatment effect (ATE) or any other averaging technique described previously with respect to Figures 2 and 8. A decision about whether or not to action the proposed treatment of the (to-be) intervened-on variable or variables x_j may be made in dependence on the estimated treatment effect. If the treatment is determined to be effective and to have a positive effect on the target variable xi (e.g. having above a threshold effect), then the treatment may be actioned; but otherwise the treatment may not be actioned. EXTENSION TO ACCOMMODATE TIME SERIES DATA Figure 12 illustrates the concept of a temporal causal graph by way of example. The modelled scenario may be divided into a series of two or more time steps, including a current time step t representing a present time, and one or more preceding time steps t-1, t-2, … representing one or more points in time in the past. Conceptually the nodes of the graph may be thought of as arranged into layers, each layer corresponding to a given time step. Each node in the graph represents a particular one of the variables of the input feature vector at a given time step. In the layer representing the present time t, the graph includes nodes representing a plurality of variables x1,t … xD,t of the feature vector at the present time (up to all the features of the feature vector). In the layer representing the first preceding time step t-1, the graph includes nodes representing a plurality of variables x_1,t-1 … x_D,t-1 at the previous time step. These may include some or all of the same variables x1,t … xD,t as included at the first time step, i.e. representing the same features of the feature vector, but instantiated at the preceding time step (the values of variables may change over time). Preferably it is the same set of variables as in the first time step. Similarly each further preceding layer before that (t-2 … t-T), if present, represents variables of the feature vector at earlier times going backward in time. Each layer may represent some or all of the same set of variables as included in other layers, but at different steps in time. In the case of more than two time steps, these may or may not be evenly spaced in time. Preferably the time steps are evenly spaced. As well as edges 1202 between nodes in the layer representing variables x1,t … xD,t at the present time t, the temporal graph is able to include one or more edges 1204 directed from the nodes representing one or more variables x_1,t-1…x_D,t-1; … x_1,t-T … x_D,t-T in layers representing one or more preceding steps in time t-1 … t-T to nodes in the layer representing variables x1,t … xD,t at the present time t. These could include an edge from a given variable at a preceding time step to the same or a different variable at the present time step. The edges could include an edge from a given variable in a preceding time step to a single other variable in the current time step, or from a given variable in a preceding time step to multiple variables at the present time step. The edges could include edges from a plurality of different variables in a given preceding time step, each to a respective set of one or more variables in the present time step, where the sets may or may not include one, some or all of the same variables. The edges could include edges from variables in a plurality of different preceding time steps, each to a respective set of one or more variables in the present time step, where again the sets may or may not include one, some or all of the same variables. In general there may be any combination of edges between any one or more variables of the feature vector at the present time step or any one or more preceding time steps and any one or more variables of the feature vector at the present time step; excepting only that edges between time steps must be directed forward in time, i.e. from a preceding time steps to the present time step. Preferably however, the graph formed by edges from present time steps to present time steps should not contain any directed cycles. (i.e. it should form a directed acyclic graph). A particular instance of a temporal graph G may be sampled from a temporal graph distribution qφ where qφ now represents a temporal graph distribution rather than a static graph distribution as specified earlier. In other words qφ now represents not only the probabilities of causal edges existing between various possible combinations of variables x_1,t … x_D,t at the present time t; but also from each of one or more of variables x1,t-1…xD,t-1; … x1,t-T … xD,t-T at one or more preceding time steps t-1…t-T each to a respective one, more or all of the variables x1,t … xD,t at the present time step t. Preferably q_φ includes an entry representing the probability of an edge existing between every possible combination variables in preceding and present time-step layers. As a temporal graph distribution, qφ may now be expressed as a matrix: qφ = {q^φt, q^φt-1, … q^φt-T}; where q^φt is the same as the static version of qφ that was defined earlier in relation to Figure 2, i.e.: ^^ _^^ _{^^ = (}

Which could also be written as q^φt = {φ_exists(i,j); φ_dir(i,j)} counting through all modelled possible combinations of i and j. The parameter φ_exists(i,j) represents the probability that a causal edge exists between x_i and x_j; and φ_dir(i,j) represents the probability that the direction of the possible edge between xi and xj is directed from xi to xj (or vice versa). Other equivalent mathematical representations of the same functionality may also be possible. Note this submatrix q^φ _t is just for present time to present time, and so does not included existence of edges {1,1}, {2,2} … {D, D}. The other submatrices q^φt-1 … q^φt-T modelling the distribution of edges from historical time steps do not contain entries for direction as edges from the past, if they exist, must always be directed forward in time. However these submatrices may include edges {1,1}, {2,2} … {D, D}. Therefore

where τ = 1… T. In other words: q^φt-τ = {φτ_exists(i,j)} counting through all modelled possible combinations of i and j. The parameter φ_τ_exists(i,j) represents the probability that an edge exists between xi at past time τ and xj at the present time t. Again other equivalent mathematical representations of the same functionality may also be possible. Figure 13 shows schematically an extension of the model of Figure 2 to accommodate the modelling of temporal causality using a temporal graph sampled from a temporal graph distribution. Based on this the extended model includes embeddings for different time steps. The extended model 104’’ is the same as the base model 104 or 104’ described in relation to Figure 2 or 10, except where specified otherwise where the extended model 104’’ builds on the base model.As well as a respective first encoder g^ei,t for each of the variables xi,t at present time t, as in the model of Figure 2, the extended model 104’’ also includes a respective first encoder g^e _i,t- _τ for each of a plurality of the variables x_i,t-τ of the feature vector for each of the one or more preceding time steps t-τ, where τ=1…T. Each first encoder g^ei,t; g^ei,t-τ is arranged to receive an input value of its respective variable xi,t ; xi,t-τ at its respective time step t, and to encode into a respective embedding e_j,t; e_{i,t- τ} for its respective combination of variable index i and time step t. The selector ^^ receives a selection of a particular selected variable of index i from among the variables (i.e. features) of the input feature vector {x₁…x_D}. The selector ^^ also takes as an input a particular sampled graph G from the graph distribution q_φ, remembering that q is now a temporal graph distribution and so G may now be a temporal graph containing one or more edges from one or more variables at a past time t- τ to one or more variables at the present time t (depending on what edges get sampled from the distribution q_φ on any given graph sampling event). The edges present in G are determined to exist or not pseudorandomly according to the corresponding probabilities specified in the temporal graph distribution qφ. The selector ^^ then selects the parents Pai that are causes of the selected variable xi in the sampled temporal graph G, including any parents from past time steps t-1…t-T. For these parents, the respective embeddings e_j,t; e_j,t-τ are generated from the respective encoder g^e, and the selector inputs these into the combiner 202 (e.g. adder or concatenator) to be combined (e.g. summed or concatenated) into the combined embedding eC for the selection variable xi , which may also be written {e_c}_I to indicate that it will vary with i. The combined embedding e_C is input into the respective decoder g^d _i associated with the selected variable in order to thereby generate a respective reconstructed value xi’, in the same way as the decoder side of Figure 2, except that this is now specifically the reconstructed value of the selected variable xi for in the present time step t. Using a temporal graph can lead to more accurate predictions of treatment effects or the like since the model 104’’ has the capacity to more accurately reflect the causation of the modelled scenario, which may be temporal in nature. In embodiments, the model 104’’ may also be used to predict the best order in which to apply treatments, the best timing in which to apply one or more treatments, and/or to predict the likely timing of an effect of one or more treatments. In example use cases, the model 104’’ (once trained) may be used to determine treatments to apply to a data centre or cloud server system, by predicting the effects of one or more server units failing or being taken offline, or having a high load. In another example, the model 104’’ may be used to determine the order in which to apply controls in an autonomous vehicle or robot, or to apply certain actions in an industrial process or farming (e.g. when to apply pesticides and fertilizer, etc.). In embodiments the reconstructed value may be combined (e.g. summed) with a noise term to generate the simulated value Xi. Alternatively the noiseless reconstructed value xi’ could simply be used as the simulated value without including noise. If noise is included, this may be combined with the reconstructed value xi’ as an additive noise term, or alternatively other non-additive means of combining a noise term with the reconstructed value xi’ may be used. Either way, in some embodiments a static noise term z_i as discussed previously could be used. However more preferably, the static noise term is replaced with a history-dependent noise term ε_i,t. The history-dependent noise term εi,t is a noise term that takes into account the values of the preceding parent variables Pai^G(<t) of the selected variable xi in the sampled graph G. This means the model will take into account the past smoothness or noisiness of time series data in which at least one variable has different values over time. Thus by including the optional history-dependent noise term, the model 104’’ will be a better specified representation of the real-world scenario being modelled, e.g. taking into account the level of fluctuation in a symptom of a patient, or a sensor reading measuring the state of a device or its environment, or traffic conditions experienced over a network, or a condition of a crop of plants, or so forth. A simulated value Xi,t that is determined based on a temporal causal graph G, and which includes history dependent noise, could be expressed as: Xi,t = fi,t [ Pai^G(< t), Pai^G(t), ^t,i ] (1) In embodiments disclosed herein, this formula is broken down and implemented as follows. X_i,t = f_i[ Pa_i ^G(< t), Pa_i ^G(t)] + g_i[ Pa_i ^G(< t), ^_t,i) ] (1a) Equation 1 is not very practical in terms of modelling as it is too general. On the other hand Equation 1a (also called equation 5 later in the subsection “Example Implementation: Auto- Regressive DECI) breaks the representation down into a part that is based only on the parents from the present and preceding time steps, but not the noise; and a part that is based only on the parents from previous time steps and the noise, but not the parents from the present time step. If the model is too general as in Equation 1, then it may not work well as it cannot be identified from observational data alone. On the other hand if the model is made too specific, then it may also not work well because sometimes the data may violate the assumptions of the equation. The inventors have found that Equation 1a provides a good workable trade-off between specificity and generality. Equation 1a may be implemented as shown in Figures 13 and 14. The embeddings ei,t-τ of the identified parents Pai^G(<t) from any preceding times steps – but not the parents Pai^G(t) from the current time step – are selected by a selector 1401 and input into a noise model 1302. In embodiments these are the same embeddings that you use to generate the reconstructed values x_i’ in Figure 13, generated by the same encoders g^ei. Alternatively however the embeddings used here in the noise model 1302 of Figure 14 could be separate embeddings generated by separate encoder networks than used in Figure 13. This selector 1401 may be implemented by some or all of the same code as the selector ^^ (it will be appreciated that the figures are somewhat schematized and do not necessarily imply that the illustrated functions have to be implemented by separate modules of code). These embeddings ei,t-τ are input into the history dependent noise model 1302 which is configured to generate the history-dependent noise term ε_i,t in dependence on the input embeddings. An example implementation of the history-dependent noise model 1302 is shown in Figure 14. The noise model comprises a combiner 1402 (e.g. adder or concatenator), a demultiplexer 1404, and a respective decoder g^d _i associated with each of the possible variables x_i…x_D that may be selected as the selected variable xi. The combiner 1402 of the noise model 1302 may share some or all of the same code with the combiner 202 of the wider model 104/104’’ as described previously with respect to Figures 2 and 13. The demultiplexer 1404 of the noise model 1302 may share some or all of the same code with the demultiplexer 204 of the wider model 104/104’’ as described previously with respect to Figures 2 and 13. In embodiments, the decoders g^di used here in the noise model 1302 of Figure 14 may be different than those used to decode the combined embedding e_C in the wider model 104/104’/104’’ as described previously with respect to Figures 2 and 13. For example the x'i decoder would only output a single value: x'i, but in embodiments the output ηi that is the input to the noise generator 1406 actually has a slightly different form than the output of the decoder that produces x'_i. E.g. Specifically, η_i might be the mean and variance (or standard deviation) of a distribution such as a Gaussian distribution (two values) and the noise generator 1406 would generate a sample from this distribution. Alternatively however, in other embodiments the same decoders g^d _i may be used both to decode η_i in the noise model 1402 of Figure 14 and to decode C_i in the wider model. Again it will be appreciated that the figures are somewhat schematized and do not necessarily imply that the components 1402, 1404 and g^di are separate modules of code than those shown in Figure 13). The noise model 1302 further comprises a noise generator 1406. The combiner 1402 combines (e.g. sums or concatenates) the embeddings ei,t-τ of the selected past parents Pai^G(<t) into a further combined embedding Ci. The demultiplexer 1404 routes this further embedding Ci into the respective decoder g^d _i associated with the selected variable x_i. This decodes the embedding C_i into a set of one or more parameter values ηi, which are the values of one or more parameters of a probabilistic distribution. The noise generator 1406 then generates the history dependent noise term ε_i,t based on the probabilistic distribution as parameterized by the one or more parameter values. E.g. the probabilistic distribution could be a Gaussian or a spline function. Figure 4 shows an example of a probabilistic distribution such as a Gaussian that may be parameterized by two parameters, a centre point (e.g. mean μ) and a width (e.g. standard distribution σ or variance σ²). In one example implementation, the noise generator 1406 could generate the distribution parameterized by the parameter(s) ηi (the mean μ and standard deviation σ), and then sample the value of the history dependent noise term εi,t directly from the generated probabilistic distribution. Alternatively however, more preferably the noise generator 1406 implements the sampling by taking a sample Ei from a fixed distribution, e.g. a unit distribution such as a unit Gaussian (with mean of 0 and standard deviation of 1), and then apply a transform function F which transforms the combination of the sample E_i from the fixed distribution and the generated parameter(s) η_i into a sample ε_i,t from the parameterized distribution (without having to actually generate the entire parameterized distribution). In other words εi,t = F(Ei, ηi). The various encoders g^ei,t; g^ei,t-τ ; decoders g^di and noise model 1302 may be implemented as modules of software as part of the model 104’’, stored in memory and arranged to run on one or more processors in a similar manner as described previously in relation to Figure 2. The extended model 104’’ may be trained in the same way as described previously with reference to Figures 2 and 7, but using a temporal graph sampled from the temporal graph distribution qφ = {q^φt, q^φt-1, … q^φt-T} instead of a static graph sampled from the static graph distribution, and including the embeddings of the past parents Pa_i ^G(<t) as well as the present parents Pa_i ^G(t), and optionally also including the history dependent noise εi,t. Figure 15 shows the method of training the extended, temporal model 104’’. At step V10 the method receives a new temporal training data point comprising values of some or all of the variables (features) of the feature vector {xi…xD} at two or more time steps. In embodiments it is all of the variables, but alternatively the model could be extended to handle missing data.At step V20 the method samples a temporal causal graph G from the temporal graph distribution q_φ = {q^φ _t, q^φ _t-1, … q^φ _t-T}. At step V30 the method selects a selected variable x_i with index i from amongst the observed variables of the input feature vector. At step V40 the method selects the instantaneous parents Pai^G(t), i.e. from the present time step t, of the selected variable x_i. At step V45 the method selects the lagged parents Pa_i ^G(<t), i.e. from the one or more preceding time steps t-1…t-T, of the selected variable xi. At step V50 the method generates the noiseless reconstructed value of xi’ of the selected variable xi, by combining the embeddings e_i,t; e_i,t-1; … e_i,t-T of the selected instantaneous and lagged parents. Optionally at step V55 the method also generates the history dependent noise εi,t. In embodiments however, only the noiseless reconstruction x’i is used during training, whereas the simulation Xi with noise is used in treatment effect estimation. At step V60 the method determines the difference between the simulated value X_i of the selected variable and the actual input value x_i. The simulated value may include the noise term, or in less preferred embodiments the noiseless reconstruction xi’ could be used as the simulated value in step V60. At step V70 the method loops back to step V30 via step V80 where it selects a new one of the observed variables of the feature vector, and repeats the method until reconstructed xi’ values have been generated for the whole feature vector. At step V90 the method then updates the parameters (e.g. weights) of all the constituent models (e.g. NNs) which include at least the encoders g^e, decoders g^d and temporal graph distribution qφ. The update is done based on the evaluating a measure of overall loss between the simulated value Xi or reconstructed value x’i of the feature vector and the input value x_i. E.g. this may be done using an ELBO function. In embodiments only the noiseless reconstruction is used during training, whereas the simulation with noise is used in treatment effect estimation. At step V100 the method loops back to V10 and repeats with a new training data point, tuning the parameters of the model 104’’ each time. This continues until the last training data point in the batch has been used, then the method proceeds to step V110 where the model 104’’ is made available to be used for making predictions, e.g. treatment effect estimation. The method of using the extended temporal model 104’’ to perform treatment effect estimation is analogous to that described earlier with reference to Figures 2 and 8. Figure 16 shows a particular example method of using the model 104’’ to estimate that treatment effect that will be experienced at a future point in time. At step W10 the method selects a target variable xA as the selected variable and a future time t’ as a target time. A here is the value of the index i for the target variable. The method also selects at least one variable xB to be intervened on (i.e. treated). This may be referred to herein as the intervened-on variable, meaning the variable for a which an intervention (treatment) is proposed, i.e. an effect of intervening on that variable is to be estimated. B is the value of the index i of the intervened on variable x_B. The method also selects a value I of the proposed intervention on x_B, and a time t~ as the time of the proposed intervention. The target time t’ is the time for which the effect on the target variable is to be estimated, i.e. the time of interest. The time t~ on the other hand is the time at which it is proposed to perform the intervention I on the intervened-on variable xB. The method begins with t=t~ where t~ is the start time, and at W20 samples a temporal graph G. At step W30 the method extracts instantaneous observations of the feature vector {x_i}, where i counts through the indices 1…D of some or all of the features of the feature vector. This may include at least those that are observed and/or treated, and/or at least those that are immediate parents of the target variable i=A in the sampled graph. For any treated (i.e. intervened-on) variables x_B these are set to their fixed, treated values I. For any observed values these are set to their observed values. In embodiments it may be assumed that observed or treated values are available for the entire feature vector at the start time t=t~. At step W35 the method extracts any past observations of the feature vector {x₁…x_D}, i.e. for times t-T…t-1. At step W50 the method uses the values of the present and previous parents to generate the reconstructed value x’i,t, cycling through the indices i=1…D of some or all of the variables as the selected variable so as to reconstruct some or all of the feature vector in the manner discussed previously. The reason for reconstructing multiple variables, even though only one target variable xA may be ultimately of interest, is that when iterated over multiple time steps then these may have a knock on effect on the target variable x_A even if not immediate parents of x_A. At optional step W40 the method generates the noise distribution parameter(s) η_i , for some or all of i=1...D, based on the historical data and sampled past parents, in the manner described previously; and then step W55 generates the history dependent noise based on this.At step W60 the method generates the simulated values {Xi,t} for i=1…D. Note: the representation of steps W20-W60 in Figure 16 may be somewhat of a simplification compared to the actual implementation. In embodiments these steps may be arranged into a sub- loop to count through i=1…D, reconstructing each variable in turn, where the order of the indices from 1 to D preferably represent a topological order of the nodes in the sampled graph. That is, in the first cycle of the subloop the value of the variable i=1 with the highest topological order is determined, and then in the next cycle of the subloop the value of variable i=2 with the next highest topological order is determined, which may depend on the reconstructed or simulated value of the first variable, and so forth. Operating in this manner means that the method can takes into account that some variables may have a knock-on effect on others (not necessarily just their immediate parents). For example at a given t, x2 might have knock-on effects on x3 (or similar), and as such the sample X_i,t can have an impact on X_j,t. The sub-loop will take account of this. This is shown in Figure 17. At step W70 the method loops back, via step W75 where t is incremented by one step, to the point of branching to steps W30/W35/W40. Each time the reconstructed or simulated values from the preceding step are fed back to the input of the model as the input values, unless intervened-on or observed values of the variables are to be used for any intermediate time steps between start time t=t~ and target time t=t’ in which case those values are used for those variables. The method continues repeating through to step W60 until future time t’ is reached. After this the method proceeds to step W80 where, if not finished, it loops back to step W20 and repeats one or more times, resampling the graph at step W20 each time. Once the method has completed all the desired graph sampling loops, it finally proceeds to step W90 are an average treatment effect for the target variable x_i of index i is determined averaged over all the sampled graphs. The temporal graph is only sampled once throughout the entire time t to t'. Since there is additional time information, both the variable index i and a future time t' are specified (t'>=t) for the target variable.Therefore, in order to simulate a variable in the future, one needs to start with the current time t and simulate the all variables at the current time t (W30 - W60). Then, one can proceed to simulate the variables at t+1. Repeat the above procedure until reaching t' (the loop of W70). This will generate simulation values of the target variable w.r.t the current temporal graph. Since it is desired to average multiple graphs, the procedure is then repeated one or more times (the loop of W80).Note also again that the flow charts may be somewhat schematized and in practice one or more steps and/or loops shown sequentially may in fact be performed in parallel. Figure 17 shows explicitly the sub-loop mentioned previously with respect to steps W20-W60 of Figure 16. Figure 17 may be considered as a more detailed expansion of Figure 16. At step W10 the method selects the variable xA with index i=A as the target variable, and selects at least one variable xB (index i=B) as the intervened-on variable, and I as the value of the intervention. The method also selects time t~ as the time of the intervention and t’ as the target future time. For example if the model 104’’ is being used to predict the effect of a heating control loop, x_A could be the temperature of a target environment or device, xB could be the thermostat setting, I the value to which the thermostat is proposed to be set (the proposed treatment), time t~ is the time of applying the thermostat setting, and time t’ is the time for which it is desired to predict the resulting change in the temperature xA. The method begins with the current time t set to an initial value, e.g. t=1, which could be the treatment time t~ or a time earlier than the treatment time. At step W20 the method samples a temporal graph G from the graph distribution q_φ. At step W25 the method sets i=1. Preferably the index i counts in topological order of the nodes in the sampled graph, such that the variable x1 indexed by i=1 has no parents, then the variable x2 indexed by i=2 has either no parents or only x₁ as a parent, and so forth. At step W27 the method determines whether the current time t equals the treatment time t~ AND the current value of the variable index i equals B, that of the intervened-on variable (i.e. treated variable) xB. If both of these conditions are met then the method branches to step W30a where it sets the treated variable xB to the value of the intervention I. Though only a single treatment index value B is shown, more generally B could represent a set of one or more index values and I could represent a set of one or more corresponding intervention values, such that one or more variables i will be set to their corresponding intervened-on values I as the method counts though i via the loop W30b. Also, though not shown, optionally step W27 could check whether there exists an observed value of the current variable i for the current time t, and if so branch to W30a where that variable is set to its observed value. In other words, more generally step W27 could represent determining whether the current variable i at the present time t has a pre-specified value determined by some other means than reconstruction by the model 104’’ (whether because it is intervened-on or observed) and step W30a may represent setting any such variable to its specified value. In some cases the method could begin with a full set of observed and/or treated variables at the initial time t=1, such that at the initial time step all the variables of the feature vector are set to specified values. If however it is determined at step W27 that the value of the present variable xi,t at index i at the present time t is not specified, then the method branches to step W29 where it extracts the values of any instantaneous parent(s) of xi,t from the sampled graph, and step W35a where it extracts the values of any lagged parents of xi,t according to the sampled graph. The method then proceeds to step W50 where it generates a noiseless reconstructed value x’_i,t of the present variable x_i,t for the present time t, using means as discussed previously in the description of the model 104’’. Optionally steps W45 and W55 are also performed to generate the history dependent noise term εi,t as also discussed previously.In this case then at step W60 the noise is added to (or otherwise combed with) the reconstructed variable x’i,t to generate the noisy simulated sample Xi,t corresponding to the present variable x_i,t. The method then proceeds to step W65 where it determines whether the counting index i has reached the value D of the index of the last variable in the feature vector (preferably sorted in topological order). If not the method branches to step W30b where it increments i by 1 and then loops back to step W25 and repeats until it has reconstructed the entire vector (or at least as much of it as is to be reconstructed). The steps W30a and W35a and the loop via W30b may be considered to correspond to the steps W30 and W35 in Figure 16, in that steps W30a sets pre- specified (treated or observed) variables to their specified values and step W35 extracts the historical values, while the loop via W30b ensures that a reconstructed value x’<i,t is generated for any unspecified variable that comes before the present variable xi,t in the counting order (e.g. topological order). This is relevant since these variables x’_<i,t could have an effect on the present variable. Counting through all variables i=1…D will mean that, over the entire count (i.e. all cycles of the subloop W27 to W30b), the past and present values of all the variables in the feature vector will be extracted as illustrated at a higher level of abstraction in steps W30 and W35 in Figure 16. If the last index D in the feature vector has been reached for the present time t, then at step W the method checks whether the present time t has reached the target time t’ yet. If not, the method branches to step W75 where t is incremented by 1, and then loops back to step W25 where it repeats with the new value of t. When the method does reach the target time t’, then the method proceeds to step W80 where it determines whether any more graphs are to be sampled. If so, the method loops back to step W20 where it repeats with a newly sampled temporal graph. Once the method has been repeated with the desired number of rounds of graph sampling, then it proceeds to step W90 where the average treatment effect on the target variable xA is determined, by averaging over the multiple sampled graphs as discussed previously. EXAMPLE USE CASES The variables of the feature vector xi = [x1…xD] may represent properties of a tangible real-world entity being modelled, or an environment in which the entity resides. The entity may comprise a living being; or a physical object such as an electronic, electrical or mechanical device or system. Alternatively or additionally the object may comprise a piece of software running or to be run on a computer system. In one example use case, the modelled real-world entity may comprise a living being, e.g. a human or animal. The disclosed model 104, 104’ or 104’’ may be used to determine whether to action a treatment to the being in order to affect a condition of the being (e.g. to cure or alleviate a disease by which the being is afflicted). In this case the target variable xi represents a measure of the condition, e.g. a symptom experienced by the living being as a result of the condition, and one or more intervened-on variables x_j represent one or more properties of the being or its environment which are susceptible to being controlled as a treatment (e.g. dosage of a drug, or a change in habit or environment, etc.). E.g. this may be used in a medical setting to treat a patient, where the potential treatment in question could be, for instance, a drug, surgery, or a lifestyle change. In another scenario the living being may comprise a plant or crop of plants, such as an agricultural crop, and the potential treatment could be the timing of sowing seeds, when to water or cover the plant/crop, or whether or when to apply a chemical such as a fertilizer or pesticide and/or in what quantity, etc. Whatever form the living being takes, a decision as to whether to perform the proposed treatment on the intervened-on variable(s) as a means to treat the target condition may be made in dependence on whether a confounder is estimated to exist between the proposed treatment and the target condition, and/or an estimate strength of the confounder. And/or, the decision may be made in dependence on the estimated treatment effect, e.g. ATE. In another example use case, the real-world entity may comprise a mechanical, electrical or electronic system or device; and the disclosed model 104, 104’ or 104’’ may be used to determine whether to action a treatment to affect a state of the system or device, such as to repair, maintain or debug the system or device. In this case the target variable xi may represent a measure of the state of the system or device, e.g. sensor data measuring wear, operating temperature, operating voltage, output or throughput, etc. The one or more intervened-on variables x_j may represent one or more properties of the device or its environment susceptible to being treated, e.g. vibration, humidity, ambient temperature, applied current or voltage level, etc. E.g. this could be used in an industrial setting to service industrial equipment such as factory equipment. A decision as to whether to perform the proposed treatment on the intervened-on variable(s) as a means to treat the target state may be made in dependence on whether a confounder is estimated to exist between the proposed treatment and the target state, and/or an estimate strength of the confounder. And/or, the decision may be made in dependence on the estimated treatment effect, e.g. ATE. In another example use case, the real-world entity being modelled may comprise software that is run, or to be run, on one or more processors in one or more computer devices at one or more locations; and the disclosed model 104, 104’ or 104’’ may be used to determine whether to action a treatment to try to optimize the running of the software. In this case the target variable xi may comprise a measure of a current state of the software, e.g. memory or processing resource usage, or latency, etc. The one or more intervened-on variables x_j may represent any property capable of affecting the running of the software, e.g. input data, or rebalancing or the load across a different combination of processors or devices. For example this could be used in a data centre or the cloud for load balancing of software run across different server devices. A decision as to whether to perform the proposed treatment on the intervened-on variable(s) as a means to optimize the running of the software may be made in dependence on whether a confounder is estimated to exist between the proposed treatment and the target state, and/or an estimate strength of the confounder. And/or, the decision may be made in dependence on the estimated treatment effect, e.g. ATE. In another example use case the real-world entity being modelled may comprise a network, e.g. a mobile cellular network, a private intranet, or part of the internet (e.g. an overlay network such as a VoIP network overlaid on the network). The disclosed model 104, 104’ or 104’’ may be used to determine whether to action a treatment to try to optimize the operation of the network. In this case the target variable x_i may represent any state of the network that it may be wished to improve, e.g. a property of network traffic such as end-to-end delay, jitter, packet loss, error rate, etc. The one or more intervened-on variables xj may represent and property capable of affecting the target variable, e.g. balancing of traffic across the network, routing or timing of traffic, encoding scheme used, etc. A decision as to whether to perform the proposed treatment on the intervened-on variable(s) as a means to optimize the network performance may be made in dependence on whether a confounder is estimated to exist between the proposed treatment and the target state, and/or an estimate strength of the confounder. And/or, the decision may be made in dependence on the estimated treatment effect, e.g. ATE. The decision(s) may be made adaptively (i.e. dynamically in response to changing conditions), or as part of network planning. In another example use case, the real-world entity may comprise an autonomous or semi- autonomous self-locomotive device such as a self-driving vehicle or robot. The disclosed model 104, 104’ or 104’’ may be used to determine whether to action a treatment in the form of a control signal to control the motion of the device. In this case the target variable x_i may represent sensor data providing information on the device’s environment, e.g. image data or other sensor data (such as distance or motion sensor data) providing information on the environment of the device (such as presence of obstacles, location of another object to interact with). The one or more intervened- on variables x_j may comprise a control signal that can be applied to control the device, e.g. to steer the vehicle or apply brakes or lights, or to move a robot arm in a certain way. A decision as to whether to perform a certain control operation (a type of “treatment”) represented by the intervened-on variable, as a means to achieve a certain effect on the environment (e.g. avoid an obstacle or interact with another object), may be made in dependence on whether a confounder is estimated to exist between the proposed variable to be controlled and the target outcome, and/or an estimate strength of the confounder. And/or, the decision may be made in dependence on the estimated treatment effect, e.g. ATE. EXAMPLE IMPLEMENTATION: AUTO-REGRESSIVE DECI Further details of some example implementations of the various concepts discussed above are now described, by way of example only. Identifying the causal relations between different variables is a question in many fields of scientific enquiry. Randomized control trial is the gold standard for discovering such relationships. However, it may not be always feasible due to practical or ethical reasons. Therefore, with no access to the interventional data, causal discovery may be used to discover the causal relations between different variables using the observational data. Recently, there has been a increasing trend towards causal discovery for time-series data. This is a relevant task in many fields of research including engineering, biology and medicine. In a temporal system, one key underlying assumption that greatly simplifies the analysis is that the arrow cannot go against the flow of time. Causal discovery methods for temporal data may be based on the framework of Granger causality, stating that a variable in the past Granger causes another in the current time if the lagged one contains some unique information about the current one. These methods can be roughly categorized into (1) vector autoregressive based methods and (2) deep learning based approaches. However, one hidden assumption of Granger causality is that the time sampling interval is frequent enough such that no aggregation effect exists. In practice, this assumption can be easily violated, where the aggregation manifest itself as the instantaneous effect, i.e. the arrow between variables can exist within a time frame. There is another line of research based on structure equation models (SEM), which can simultaneously model the lagged and instantaneous effect. The linear vector auto-regressive may be combined with non-Gaussian noise for causal discovery. For example a general structure equation model for time-series data, called TiMINo, provides a general theoretical framework for structure identifiability. The above methods assume that the noises are mutually independent and stationary in a way such that their distributions have independent (and learnable) parameters. However, this may not hold in many scenarios. For example, under education context, the observational noise of whether a student correctly answers a question should depend on his/her past learning history. If the student correctly answered similar past questions in a consistent way, the observational noise should be small. On the other hand, if the answer history is similar to a random guess, the observational noise should be large. This history-dependent noise distribution is ubiquitous in real life and cannot modelled by the aforementioned SEM. To address the aforementioned issue, a particular embodiment of the present disclosure provides a novel SEM for time-series data based on the framework known as deep end-to-end causal inference (DECI) (Geffner et al, 2022). This may be named herein auto-regressive DECI (AR- DECI), which can model both non-linear lagged and instantaneous effect with history-dependent noise distribution. To implement this, the following contributions are provided. ^ A new formulation of temporal SEM with history-dependent noise, called AR-DECI, which combines SEM with vector auto-regressive to model both non-linear lagged and instantaneous effects. AR-DECI also adopts a Bayesian view of graph learning by using variational inference to approximate the graph posterior, which provides uncertainties over graphs under limited data. Also, we show one can compute interested treatment effect estimation by leveraging the fitted AR-DECI, e.g. conditional averaged treatment effect (CATE) with time- invariant interventions. ^ Theoretically, we show that the proposed AR-DECI is structure identifiable under assumptions. To achieve this, we provide a general framework for showing structure identifiability with history-dependent noise, where AR-DECI is a special case. Furthermore, we show AR-DECI unifies several aforementioned approaches based on SEM and Granger causality. First we will briefly introduce necessary preliminaries required for building AR-DECI model. In particular, we will focus on the basic concepts of structure equation model, DECI (Geffner et al. 2022), and typical approaches for causal time series: Granger causality and vector auto-regressive models. Note that capital X may now be used to represent both input values and simulated values of the variables. Structure Equation Models (SEM): Consider a multivariate time-series X_t = (Xⁱ _t)_i∈_V where V is a set of nodes and |V| = D. SEM describes the functional relationships between individual variable Xⁱt across and within the time frame given a graph G. In particular, a general form of SEM is

where Pa_i ^G(< t) contains the parent node specified by G in previous time frame (lagged parents) (i.e. parent nodes at time t−1, t−2, ...), Pa_i ^G(t) is the parent node at current time t (instantaneous parents), ^^t,i is the exogenous noise that is mutually independent to other variables in the model and fi,t describes the functional relationships between the parents and child node Xⁱt . The above SEM induces a joint probability distribution over the entire time series {X_t}_t=0…T. Throughout the following, we make two underlying assumptions for SEM: (1) the function f_i for all i ∈ V and t is differentiable w.r.t all the arguments and (2) the induced joint distribution is absolutely continuous w.r.t the Lebesgue measure. Although the above form is very general, it is not very practical to be directly used for causal discovery due to the non-identifiability issue. Deep End-to-end Causal Inference (DECI): Recently, Geffner et al. (2022) disclosed an end- to-end framework, called DECI, for causal inference with static data. DECI is based on a specific type of SEM, name additive noise model (ANM), which is structure identifiable under mild conditions. Specifically, the SEM of DECI is defined as: Xi = fi [ PaG(Xi) ] + ^i (2) where we drop its dependence of time t compared to Eq. (1) due to its static nature, and ^i is the additive independence noise. DECI also adopts a Bayesian view for graph learning by applying variational inference to approximate the graph posterior p(G|D) where D is the training dataset. By leveraging the trained DECI, one can not only perform causal discovery but also inference tasks such as estimating (conditional) average treatment effect, which provides us an end-to-end pipeline from observation data to the interested causal quantities. Granger Causality: Granger causality has been extensively used as the underlying principle for identifying causal relationships in time series data. It is based on the idea that Xⁱ does not Granger causes X^jt if the history of Xⁱ does not help the prediction of X^jt given the past of all other time series Xk for k ≠ i. Formally, we can express Granger causality under the SEM framework: Definition 1 (Granger Causality): Given a multivariate stationary time series {X_t}_t=0…T and a SEM fj defined as: X^jt = fj,t [ PaⁱG(< t), ^t,j ] (3) where Paⁱ _G(< t) is the lagged parents of node X^j _t and ^_t,j is the independent noise. Under this setup, time-series i Granger causes j if ∃l ∈ [1, t] such that Xⁱ _t−l ∈ Paⁱ _G(< t) and f_j,t depends on Xⁱ _t−l. It has been shown that Granger causality is equivalent to causal relations for directed acyclic graph (DAG) if there are no latent confounders and instantaneous effect. One shortcoming of Granger causality is its incapability of handling the instantaneous effect, which can happen when the time sampling interval is not frequent enough. Vector Auto-regressive Model: To overcome the aforementioned issue of Granger causality, another line of research focuses on directly fitting the identifiable SEM to the observational data with instantaneous parents. One typical approach is called vector auto-regressive model with the SEM defined as:

where β^j is the offset, τ is the lag, B^τ ∈ ℝ^D×D is the weighted adjacency matrix specifying the connections in the DAG (i.e. if B^τ _ji = 0 means no connection between X^j _t−τ to Xⁱ _t ) and ^_t,i is the independent noise. Under certain assumptions, the above SEM has been shown to be structure identifiable, which validates its usage in causal discovery/inference tasks. However, in real life, linear relationship is too restricted to capture the complex data generating mechanism and noise distribution may also be history-depended as described earlier. Therefore, a crucial research question is how to design an SEM for such effects with structure identifiability guarantees. AR-DECI: Auto-Regressive Deep End-To-End Causal Inference In this subsection, we will first focus on formulating the SEM of AR-DECI, followed by disclosing a variational objective with continuous DAG relaxation for training. In the end, we will discuss how to estimate CATE with time-invariant interventions. Model formulation: For a multivariate stationary time-series {Xt}t=0…T with a set of nodes V and |V| = D, we assume it follows a temporal graph adjacency matrix G_0:K where K is the maximum lag K. We define G_τ,ij = 1 means there is a directed connection Xⁱ _t−τ → X^j _t and 0 means no connection. In the following, we interchange the usage of the notation G and G0:K for brevity. We provide the following SEM for AR-DECI:

where fi, gi are general differentiable non-linear functions, PaⁱG(< t) is the lagged parents of node Xⁱ _t specified by adjacency G_1:K, Paⁱ _G(t) is the corresponding instantaneous parents specified by G0, and ^t,i is the independent noise. We also assume the above model satisfies causal Markov and minimality conditions. The above SEM has a similar additive structure as ANM, which has been leveraged in previous work to develop causal time-series models. A distinction of AR-DECI is that the noise ^_t,i is transformed by a non-linear history-dependent function g_i, meaning that the transformed noise distribution is history dependent and may not be independent to the other variables. For the detailed design of the non-linear function fi, it should respect the relationships specified by graph G0:K. Namely, if X^jt−τ ∉ PaⁱG(< t) ∪ PaⁱG(t), then ∂fi/∂X^jt−τ = 0. Building on the design of DECI, we disclose:

where ζi and ℓτi (i = 1, ... ,D and τ = 0, ... ,K) are MLPs. To reduce the number of neural networks, we adopts a weight sharing scheme across nodes and lags: ζ_i(·) = ζ(·,u_0,i) and ℓ_τj(·)=ℓ(·,u_τ,j), where u_τ,i is the trainable embedding for node i at time t − τ . For the design of the noise transformation gi, it allows very flexible design choices. For the ease of evaluation and training, we limit our choice to be a conditional normalizing flow. Specifically, we choose conditional spline flow for g_i with Gaussian noise ^_t,i for all t and i. Due to the invertibility of gi, we have:

where we omit the dependency of gi with PaⁱG(< t) for clarity. Variational Inference for AR-DECI AR-DECI adopts a Bayesian view of causal discovery, which aims to learn a graph posterior distributions instead of inferring a single graph. For N observed multivariate time series X⁽¹⁾0:T , . .. , the joint likelihood of AR-DECI is:

where θ is the model parameter. Once fitted, the posterior p(G|X⁽¹⁾ _0:T , .. . ,X^(N) _0:T) incorporates the belief of the underlying causal relationships. Graph Prior: To design a proper graph prior, we aim at three criteria: (1) constrain the graph to be a DAG; (2) favour sparse graph; (3) support prior knowledge. We choose the prior as:

where h(G) = tr(e^G⊙G) − D is the DAG penalty proposed in and is 0 if and only if G is a DAG; ⊙ is the Hadamard product; G^p is the prior knowledge of the graph; λs, λp specifies the strength of the graph sparseness and prior knowledge, respectively; α, ρ characterize the strength of the DAG penalty. Since the connections from the history node to the current one can only follow the direction of the time flow, only the instantaneous connections can violate the DAG constraint. Thus, DAG penalty is only applied to G₀. Variational Objective: Unfortunately from the Bayes’ rule, the exact graph posterior p(G|X⁽¹⁾ _0:T , . . . ,X^(N)0:T) is intractable due to the large combinatorial space of all DAGs. To overcome this challenge, we adopt variational inference, which uses variational distribution qϕ(G) to approximate the true posterior. We choose independent Bernoulli distribution for each directed edge as qϕ(G). We derive the evidence lower bound (ELBO) as:

where H(q_ϕ(G)) is the entropy of the variational distribution. From the causal Markov assumption and auto-regressive nature, we can derive:

^and

(12) where z^i,(n) _t = X^i,(n) _t −f_i[ Paⁱ _G(< t), Paⁱ _G(t) ] and p_gi is defined in Eq. (7). The parameters θ, ϕ are optimized by maximizing the ELBO, where a Gumbel-softmax trick is used to estimate the gradients w.r.t ϕ. Treatment effect estimation: We now show how to leverage the fitted AR-DECI for estimating the CATE. For simplicity, we only consider a special case of CATE defined as:

We assume the conditioning variable can only be X_<t (i.e. the entire history before t), and the intervention and target variable can only be either at current time t or sometime in the future t + τ. This formulation is for simplicity, and AR-DECI can be easily extended to more general cases using the CATE tricks in DECI (Geffner et al). Once fitted, the idea is to draw target samples X^Y _t+τ from the interventional distribution p(X^Y _t+τ |X_<t, do(X^I _t),G) for each graph sample G ∼ q_ϕ(G). Then, unbiased Monte Carlo estimation can be used to compute CATE. For sampling from the interventional distribution, we can use the “mutilated” graph G_do(X^It) to replace G, where all incoming edges to X^It are removed. The intervention samples can be obtained by simulating the AR-DECI with history X_<t, X^I _t = a or b and G_do(X^I _t). EXAMPLE IMPLEMENTATION: CONFOUNDED DECI Further details of some example implementations of the various concepts discussed above are now described, by way of example only. As discussed earlier, in practice it is not always possible to directly perform interventional studies (i.e., randomised AB testings) to estimate treatment effects, and it may be required to use observational studies. That is, the causal graph and the associated causal effects are to be inferred based on observational data, in embodiments solely based on observational data. Therefore, this would appear to require an assumption that all the data provided by users already contains all the information that is needed (i.e. there should not be any unmeasured confounders). However, this assumption is often unrealistic, as the modelled scenario is also impacted by certain variables that cannot be directly measured. Therefore, it would be desirable to provide an effective and theoretically principled strategy for handling latent confoundings when performing joint causal discovery & inference. However, the existence of unmeasured confounding poses questions for causal discovery since there might exist multiple contradicting causal structures that are compatible with observations. Traditional methods for causal discovery mainly focus on estimating the causal relationships between observed variables and does not recover the causal relationships that involve unmeasured confounders. There also exists several related methods for causal discovery with latent confounders, but each comes with certain limitations on assumptions.Following the previous work by Geffner at al of developing deep end-to-end causal inference (DECI), this framework can be extended in order to handle latent confounders. The outcomes of this is to implement latent variable variants/baselines of DECI that explicitly take latent confounders into account. These variants/baselines should preferably: 1), be able to perform joint graph learning and functional learning; 2), have computational costs comparable to DECI; and 3), subsume DECI, in the sense that it should perform similarly to DECI when there is no latent confounders. To recap the causal discovery without latent confounders using DECI, recall that in DECI, we perform causal discovery under no latent confounder assumptions, but modelling the causal graph G on all observed vertexes X jointly with the observations x¹, ... , x^N as:

Here, p(G) is a prior over DAG, implemented by:

where G is the adjacency matrix representation, and ℎ⁽ ^^⁾ = ^^ ^^⁽ ^^ ^{^^ ^ ^^)} − ^^ (3) is the DAG penalty. We aim to fit θ, the parameters of our non-linear ANM, using observational data by maximizing (an lower bound of) log pθ(x¹, ... , x^N). Given θ, the question of causal discovery can be simply answered by the posterior, pθ(G|x¹, ... , x^N). In the original DECI, we maximize the ELBO of log p_θ(x¹, ... , x^N) and use mean-field approximation qϕ(G) to perform approximate inference. Causal discovery under ADMGs: in the present disclosure, we will extend DECI’s framework to consider unmeasured confounders. Unfortunately, if some variables in a DAG is unobserved, the resulting set of (conditional) independences no longer corresponds to a DAG (on the observed variables). To deal with unobserved variables, a common extension is the so called acyclic directed mixed graphs (ADMGs). In a ADMGs, a bidirected edge x ↔ y indicates the existence of a unmeasured founder u_xy, such that x ←u_xy → y. Note that this would be different than a undirected edge x−y used in (CP)DAG, which indicates either x ← y or x → y. Finally, the term acyclic in ADMG means that there is no cycles that contains only directed edges. Therefore, causal discovery under latent confounders can be formulated as the same equation as the original DECI, ^^ _^^( ^^¹, … , ^^ ^{^^} , ^^) = ^^( ^^)∏ ^^ _^^( ^^ ^{^^}| ^^) ^^ (4) except that the prior p(G) will be based on the constraint that G is a ADMG (instead of a DAG). We will discuss how to parameterize the ADMG constraint and how to parameterize p_θ(x¹, ... , x^N, G) later. Note that in principle, an ADMG causal discovery algorithm is able to express search results using (a subset of) the following structures: ^ x → y: x is the cause of y; ^ x ↔ y: indicates the existence of a unmeasured founder uxy, such that x ← uxy → y. ^ x − y: either x ← y or x → y. ^ x ⇒ y: either x ↔ y, or x → y. In the following, we mainly consider representing our posteriors using the first two, where they have better granularity than the last two. D-DECI (deconfounded DECI) Magnification-based representation: We perform causal discovery on the highest granularity, i.e. explicitly model the latent confounders (denoted by U), and learn a DAG G’ on X ∪U. In other words, G’ (on X ∪U) is a magnification of G (on X ); and the marginalized graph G over observed variables X will then be a ADMG. To relate G’ to the previous notation of G₁ and G₂, G₁ is the adjacency matrix on X; and G₂ is the adjacency matrix on U. G' is basically a big matrix on both X and U, obtained by concatenating and unpacking G1 and G2. In other words, G' = M(G1, G2), where M(.) is a pre-known mathematical function. G' may also be called magnified matrix. In short, it is obtained by concatenating and aggregating G₁ and G₂. For example, for a graph x₁-- >x2<-->x3, the corresponding matrices will be:

where the 4th row of G' corresponds to the 4th variable (additional latent confounder) implied in the bidirected edge x₂<-->x₃. That is, the 4th variable should point to both x₂ and x₃, hence those two 1s in the second and third cell of the 4th row. Under this idea, we model the causal graph G’ jointly with the observations x¹, ... , x^N, as well as latent confounders u¹, ... , u^N:

We will define constraints on G’ and the parameterization of p_θ(xⁿ, uⁿ|G’) later. We call Eq. (5) magnification-based parameterization, since the original ADMGs formulation does not explicitly assume any functional forms. To optimize θ and perform approximate inference on G’, we define: 1) A variational distribution q_ϕ(G) to approximate the intractable posterior p_θ(G|x¹, ... , x^N). Similar to the original DECI, this can be parametrized as the product of independent Bernoulli distributions for each potential directed edge in G’. We parametrize edge existence and edge orientation separately, using the ENCO parametrization. 2) An inference network q_λ(u|x) to approximate the intractable posterior p_θ, G’(u|x), e.g., a Gaussian inference network used in VAEs. The inference network qλ(u|x) has also been referred to as H earlier. We use them to build the ELBO, given by:

By committing to optimize Eq. (6), we implicitly make the following assumptions: ^ We already know the exact number of unmeasured variables. That is, the latent vertexes U is given and fixed. ^ For all unmeasured variables in U, they only possess out-going edges. That is, all variables U serve as confounders. This in practice and be implemented by model.set graph constraint() in DECI code base. Prior over ADMGs: The ADMG constraint on the graph G’ can be imposed in two ways: ^ Directly apply the DAG constraint on G’. That is,

This ensures that once G’ can be converted to an ADMG (by marginalizing U and only keeping X). OR ^ First marginalize G’ to the graph G on X , and then apply ADMG constraints on G. This can be done via the following three ADMG constraints: o Ancestral ADMG constraint:

where Gdirected,i,j = 1 if and only if xi ← xj ; and Gbidirected,i,j = 1 = G if and only if xi ↔ xj . Arid ADMG constraint: ℎ⁽ ^^⁾ = ^^ ^^⁽ ^^ ^{^^)} − ^^ + ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^⁽ ^^ _{^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^}, ^^ _{^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^ ^^} ⁾ = 0 (9) o Bow-free constraint:

By using any of those constraints, p(G’) can be implemented by:

Parameterization of pθ(xⁿ, uⁿ|G’): Similar to DECI, we can parameterize p_θ(xⁿ, uⁿ|G’) using auto-regressive ANM. That is, z = g_G’(x, u; θ) = [x, u] – f_G’ Q ([x, u]; θ). Where z are independent noise variables following the distribution Πi pzi. Then we can write the observational likelihood as:

where we omitted the Jacobian-determinant term because it is always equal to one for DAGs G’. Structural identifiability under latent confounders: Preferably, we would like to find a set of suitable assumptions such that D-DECI should satisfy the structural identifiability definition, as follows. Definition 1 (D-DECI Structural Identifiability). For a distribution pθ(x;G’) := ∫u pθ(x, u|G’)du, the graph G’ is said to be structural identifiable from pθ(x;G’) if there exists no other distribution p_θ’(x;G’’) such that G’ ≠ G’’ and p_θ(x;G’) =p_θ’ (x;G’’). Generally, such structural identifiability does not hold. However, under the constraint of bow- freeness, we will be prove identifiability, which will be discussed in the next section. D-DECI identifiability Notation and assumptions: Let x ∈ ℝ^D denote the observed data, u ∈ ℝ^P denote the unobserved data, v = [x; u]^T , and p(x, u; ^^⁰) the ground truth data generating distribution, where ^^⁰ is a binary adjacency matrix representing the true causal DAG. DDECI uses the additive noise model (ANM) defining the structural assignment x_j’ = f_j(pa(j); θ)+n_j where pa(j) denotes the parents of x_j. The function f corresponds to the decoders g^d described earlier. We’ll also use the notation pax(j) and pau(j) to denote the parents of xj in the observed and unobserved set, respectively. nj are mutually independent noise variables. The noise term n has also been referred to as z earlier. Assumption 1 (Form of fj): We assume that the functions fj are non-linear and decompose as: ^^ _^^( ^^ ^^ _^^( ^^), ^^ ^^ _^^( ^^); ^^) = ^^ _{^^, ^^}( ^^ ^^ _^^( ^^); ^^) + ^^ _{^^, ^^}( ^^ ^^ _^^( ^^); ^^) (13) Assumption 2 (Unobserved variables are latent confounders): All unobserved variables are latent confounders of observed variables, i.e. have the causal structure x_j ← u_k → x_i. Justification of this is given by the proposition that any latent variable can be reparameterised by absorbing exogenous noise to make it a latent confounder. A distribution p(x) satisfies minimality w.r.t. ^^ if it is Markov w.r.t. ^^, but no to any proper subgraph of ^^. Since we cannot access the true data generating distribution p(x, u; ^^⁰), only the marginal distribution over observed variables p(x; ^^⁰), we require a slightly different definition of minimality. We define this below. Definition 1 (Minimality): p(x; ^^) is minimal if it does not satisfy the local Markov condition with respect to the ADMG of any subgraph of ^^ . This form of minimality implies the bow-free constraint. To see this, note that the causal structure xi ← uk → xj together with xi ← xj, and xi ← uk → xj alone imply the same conditional independencies. Since the latter is a subgraph of the former, only it satisfies the minimality condition. Assumption 3 (Minimality): For a distribution p(x; ^^) generated by D-DECI, we assume the minimality condition holds. Comments on assumption of latent variables: Treatment effect is invariant to parameterisations if intervening on observed variables only. Proof of structural identifiability: Here we prove structural identifiability of ^^ given the assumptions of the previous section. Lemma 1: Let M denote a set of variables not containing v_i, and let s(v_i) denote a linear or nonlinear function of vi. Then, the residual of s(vi) regressed on M cannot be independent of ni:

Proof: Assume that [s(v_i) − g_i(M) n_i] holds, then M must contain at least one descendent of v_i as it must have dependence on the noise ni to cancel effect of ni in s(vi). We can express gi(M) as

. Since g_i operates on variables defined by non-linear transformations of the exogenous noise terms, we cannot express u_i as a_i

+ b_i(n_i). M contains a descendent of v_i, so ^^ _^^ ∪ ^^ _^^ includes at least one noise term n_k that satisfies v_i n_k (i.e. is not in v_i). Thus, terms containing ni cannot be fully removed from s(v_i) − g_i(M) and so [s(v_i) − g_i(M) n_i] does not hold. Lemma 2: If an only if is satisfied, there is a latent confounder, and no direct cause, between xi and xj . (15) Proof: Define gi and gj as gi(M) = fi,x(pax(i)) + g’i(M) and gj(N) = fj,x(pax(j)) + g’i(N) respectively. Then, (15) becomes equivalent to:

This is equivalent to: (17) And since n_i n_j, we have: (18) The first implies the existence of an unobserved mediator between x_j and x_i, the second implies the existence of an unobserved mediator between xi and xj, and the third implies the existence of an unobserved confounder. Given Assumptions 2 and 3, this indicates the presence of a latent confounder and no direct cause between x_i and x_j . Lemma 3. If (19) is satisfied, there is no causal relationship between x_i and x_j and no latent confounder.

Lemma 4. If (20) and (21) are satisfied, x_j is a direct cause of x_i and there is no latent confounder.

Thus, for any pair of variables we can determine whether there exists a latent confounder, and if not we can determine whether or not there is a direct cause and the orientation. MLE recovers ground proof: We can use the exact same proof as in the DECI paper to show that MLE recovers ground truth—we can just replace ^^ with the ADMG corresponding to ^^ projected onto x. D-DECI recovers ground truth: In the infinite data limit, we have: (22) The posterior approximation q(un|xn) should preferably be very accurate (i.e.

|| ^^⁽ ^^ _^^ ^| ^^ _^^, ^^⁾]] ≈ 0) to recover the ground truth causal structure. FINAL REMARKS It will be appreciated that the above embodiments have been described by way of example only. More generally, according to a first aspect disclosed herein there is provided (Statement 1): a computer-implemented method comprising: A) selecting a selected variable from among variables of a feature vector; B) sampling a temporal causal graph from a temporal graph distribution, the temporal graph distribution specifying probabilities of directed causal edges between different ones of the variables of the feature vector at a present time step, and from one of the variables of the feature vector at a preceding time step to one of the variables of the feature vector at the present time step; C) from among of the variables of the feature vector, identifying a present parent which is a cause of the selected variable in the present time step according to the temporal causal graph, and identifying a preceding parent which is a cause of the selected variable from the preceding time step according to the temporal causal graph; D) inputting an input value of each of the identified present and preceding parent into a respective encoder, resulting in a respective embedding of each of the present and preceding parents; E) combining the embeddings of the present and preceding parents, resulting in a combined embedding; and F) inputting the combined embedding into a decoder associated with the selected variable, resulting in a reconstructed value of the selected variable. In embodiments the method may optionally further comprise features as set out in any of the following Statements. Statement 2. The method of Statement 1, wherein: the temporal causal graph distribution sampled in B) specifies probabilities of directed causal edges existing from each of a plurality of variables of the feature vector from one or more preceding time steps, each to a variable of the feature vector at the present time step; and C) comprises, from among of the variables of the feature vector, identifying each present parent that is a cause of the selected variable in the present time step according to the temporal causal graph, and identifying each preceding parent variable which is a cause of the selected viable from a preceding time step according to the temporal causal graph. Statement 3. The method of Statement 2, further comprising: generating a history dependent noise term based on embeddings of the preceding parents; and combining the history dependent noise term with the reconstructed value of the selected variable, resulting in a simulated value of the reconstructed variable. Statement 4. The method of Statement 3, wherein embeddings of the present parents are not input into the noise model. Statement 5. The method of Statement 4, wherein the generating of the history dependent noise term comprises: combining the embeddings of the preceding parents, resulting in a further embedding; inputting the further embedding into a decoder associated with the selected variable, resulting in one or more parameter values; and generating the history dependent noise term from the probabilistic distribution as parameterized by the one or more parameter values. Statement 6. The method of any preceding Statement, wherein: B) further comprises sampling a second causal graph from a second graph distribution, the second causal graph modelling presence of possible confounders, a confounder being an unobserved cause of both of two variables in the feature vector; C) further comprises, from among of the variables of the feature vector, identifying a parent variable which is a cause of the selected variable according to the first causal graph, and which together with the selected variable forms a confounded pair having a respective confounder being a cause of both according to the second causal graph; and D) further comprises inputting the input value of the parent variable and an input value of the selected variable into an inference network, resulting in a latent value modelling the respective confounder of the confounded pair, and inputting the latent value into a second encoder, resulting in an embedding of the confounder of the confounded pair; and in E) the combining includes combining the embedding of the present and preceding parents with the embedding of the confounder of the confounded pair, thereby resulting in said combined embedding. Statement 7. The method of any of Statements 2 to 5, or Statement 6 when dependent on any of Statements 2 to 5, comprising: i) for a given training data point comprising a given combination of input values of the variables of the feature vector at the present times step and each preceding time step, repeating A)-F) over multiple selections, each selection selecting a different one of the variables of the feature vector as the selected variable thereby resulting in a respective reconstructed value, the multiple selections together thereby resulting in a reconstructed version of the training data point comprising the reconstructed values for the training data point; ii) evaluating a measure of difference between the training data point and the reconstructed version; and iii) training model parameters of the encoders, decoders, and temporal graph distributions, based on the evaluated measure. Statement 8. The method of Statement 7, comprising repeating i)-iii) over multiple input data points, each comprising a different combination of input values of the variables of the feature vector. Statement 9. The method of Statement 7 or 8, wherein: each selection in i) further comprises a history dependent noise model generating a respective history dependent noise term based on embeddings of the preceding parents; and wherein the training further comprises training model parameters the history dependent noise model based on the evaluated measure. Statement 10. The method of Statement 9, wherein in each selection in i), the generating of the history dependent noise term by the history-dependent noise model comprises: combining the embeddings of the preceding parents, resulting in a respective further embedding; inputting the respective further embedding into the decoder associated with the selected variable, resulting in a one or more parameter values; and generating the respective history dependent noise term based on the one or more parameter values. Statement 11. The method of Statement 10, comprising repeating i)-iii) over multiple input data points, each comprising a different combination of input values of the variables of the feature vector. Statement 12. The method of Statement 2, or any of Statements 3 to 11 when dependent on Statement 2, comprising: I) setting the input value of an intervened-on variable of the feature vector, other than the selected variable, to a specified value; and II) estimating an effect of the intervened-on variable on the selected variable based on the reconstructed value of the selected variable. Statement 13. The method of Statement 12, wherein: I) comprises setting the input value of a plurality of intervened-on variables of the feature vector, other than the selected variable, to respective specified values; and II) comprises estimating the effect of the plurality of intervened- on variables based on the reconstructed value of the selected variable. Statement 14. The method of Statement 13, wherein: I) comprises repeating A)-F) over multiple rounds, each using the same selected variable and setting the intervened-on variable to the same specified value, but performing the sampling of the temporal causal graph afresh each round, thereby resulting in a respective reconstructed value of the selected variable from each round; and II) comprises determining an average treatment effect of the intervened-on variable on the selected variable averaged over the multiple rounds. Statement 15. The method of Statement 14, wherein each graph-sampling event further comprises: generating a history dependent noise term based on the embeddings of the preceding parents; and combining the history dependent noise term with the reconstructed value of the selected variable, resulting in a simulated value of the reconstructed variable; wherein said estimating based on the reconstructed values in II) comprises: estimating the effect of the intervened-on variable on the selected variable based on the simulated values. Statement 16. The method of Statement 14 or 15, wherein one or both of: each round sets a plurality of intervened-on variables and/or observed variables of the feature vector to specified values (and in the case of multiple intervened-on variables, II comprises determining the average treatment effect of the plurality of intervened-on variables on the selected variable); and/or each round comprises an interior loop of D)-F) repeated around multiple iterations with the same sampled graph but for incrementing values of the present time step with each iteration, wherein each but the last iteration comprises, in addition to the target variable, reconstructing or simulating values of one or more further variables of the feature vector, and feeding back the simulated or reconstructed values as the specified input values of the next iteration. Statement 17. The method of any of Statements 12 to 16, wherein the intervened-on variable models a treatment on a real-world entity or an environment thereof and the target variable models an effect of the treatment applied to the real-world entity, and the method further comprises actioning the treatment on the real-world entity in dependence on the estimated treatment effect. Statement 18. The method of Statement 17, wherein one of: the real-world entity comprises a living being, and the treatment comprises a medical treatment to the living being or an environment thereof; or the real-world entity comprises a mechanical, electrical or electronic device or system, or an environment thereof; and the treatment comprises an act of maintaining, debugging, upgrading or controlling the device or system, or controlling the environment thereof; or the real-world entity comprises a network or software, and the treatment comprises an act of controlling the network or software. Statement 19. A system comprising: processing apparatus comprising one or more processors; and memory comprising one or more memory units, wherein the memory stores: a machine learning model comprising the encoder and decoders of any preceding Statement, and optionally the noise model and/or a inference network; and code arranged to run on the processing apparatus and being configured so as when run to perform the method of any of Statements 1 to 18. Statement 20. A computer program embodied on non-transitory computer-readable storage, wherein the computer program comprises instructions configured so as when run on one or more processors to perform the method of any of Statements 1 to 18. According to a second aspect disclosed herein there is provided (Statement 1A) a computer- implemented method comprising: a) selecting a selected variable from among variables of a feature vector; b) sampling a first causal graph from a first graph distribution and sampling a second causal graph from a second graph distribution, the first causal graph modelling causation between variables in the feature vector, and the second causal graph modelling presence of possible confounders, a confounder being an unobserved cause of both of two variables in the feature vector; c) from among of the variables of the feature vector, identifying a parent variable which is a cause of the selected variable according to the first causal graph, and which together with the selected variable forms a confounded pair having a respective confounder being a cause of both according to the second causal graph; d) inputting an input value of the parent variable into a first encoder, resulting in a respective embedding of the parent variable; e) inputting at least the input value of the parent variable (and optionally also an input value of the selected variable and/or an input value of one or more others of the variables of the feature vector) into an inference network, resulting in a latent value modelling the respective confounder of the confounded pair, and inputting the latent value into a second encoder, resulting in an embedding of the confounder of the confounded pair; f) combining the embedding of the parent variable with the embedding of the confounder of the confounded pair, resulting in a combined embedding; and g) inputting the combined embedding into a decoder, resulting in a reconstructed value of the selected variable. In embodiments there may optionally be provided a method in accordance with any of the following statements. Statement 2A. The method of Statement 1A, wherein c) comprises: identifying each of the variables of the feature vector that is a respective parent variable of the selected variable according to the first causal graph, and identifying each respective confounded pair that comprises the selected variable and a respective one of the identified parent variables according to the second causal graph; d) comprises: inputting a respective input value of each of the parent variables into a respective first encoder for the parent, resulting in a respective embedding of each respective parent variable; e) comprises: for each identified confounded pair, inputting the input value of the respective parent variable and the input value of the selected variable into the inference network, thereby resulting in a respective latent value of each of the respective confounders, and inputting each of the latent values into a respective second encoder for the respective latent value, resulting in a respective embedding of each of the respective confounders; f) comprises combining the embeddings of all the identified parent variables and confounders together, thereby resulting in said combined embedding; and g) comprises inputting the combined embedding into a respective decoder associated with the selected variable, resulting in a reconstructed value of the selected variable. Statement 3A. The method of Statement 2A, comprising: i) for a given training data point comprising a given combination of input values of the variables of the feature vector, repeating a)-g) over multiple selections, each selection selecting a different one of the variables of the feature vector as the selected variable thereby resulting in a respective reconstructed value, the multiple selections together thereby resulting in a respective reconstructed version of the training data point comprising the reconstructed values for the training data point; ii) evaluating a measure of difference between the training data point and the reconstructed version; and iii) training parameters of the inference network, first and second encoders, decoders, and first and second graph distributions, based on the evaluated measure. Statement 4A. The method of Statement 3A, comprising repeating i)-iii) over multiple input data points, each comprising a different combination of input values of the variables of the feature vector. Statement 5A. The method of Statement 4A, comprising, after the training over the multiple training data points, observing the second graph distribution to estimate whether a confounder exists between a target pair of the variables of the feature vector . Statement 6A. The method of Statement 5A, comprising, after the training over the multiple training data points, observing the latent value of the respective confounder between the target pair of variables, resulting from the inference network, as an estimated value of the respective confounder. Statement 7A. The method of Statement 5A or 6A, wherein the variables of the feature vector model properties of a real-world entity or an environment thereof, and the method comprises, in dependence in said estimating of whether a confounder exists between the target pair of variables, determining whether or not to apply a treatment to one of said target pair of variables to affect the other of said target pair of variables. Statement 8A. The method of Statement 7A, wherein one of: the real-world entity comprises a living being, and the treatment comprises a medical treatment to the living being or the environment thereof; or the real-world entity comprises a mechanical, electrical or electronic device or system, or an environment thereof; and the treatment comprises an act of maintaining, debugging, upgrading or controlling the device or system, or controlling the environment thereof; or the real-world entity comprises a network or software, and the treatment comprises an act of controlling the network or software. Statement 9A. The method of any of Statements 2A to 8A, comprising: I) setting the input value of an intervened-on variable of the feature vector, other than the selected variable, to a specified value; and II) estimating an effect of the intervened-on variable on the selected variable based on the reconstructed value of the selected variable. Statement 10A. The method of Statement 9A, wherein: I) comprises setting the input value of a plurality of intervened-on variables of the feature vector, other than the selected variable, to respective specified values; and II) comprises estimating the effect of the plurality of intervened- on variables based on the reconstructed value of the selected variable. Statement 11A. The method of Statement 9A or 10A, wherein: I) comprises repeating a)-g) over multiple graph-sampling events, each using the same selected variable and setting the intervened- on variable to the same specified value, but performing the sampling of the first and second causal graphs afresh each time, thereby resulting in a respective reconstructed value of the selected variable from each graph-sampling events; and II) comprises determining an average treatment affect of the intervened-on variable on the selected variable averaged over the multiple graph- sampling events. Statement 12A. The method of Statement 11A, wherein each graph-sampling event sets a plurality of intervened-on values to specified values, the same values each time; and II) comprises determining the average treatment effect of the plurality of intervened-on variables on the selected variable. Statement 13A. The method of Statement 9A, wherein the intervened-on variable models a treatment on a real-world entity or an environment thereof and the target variable models an effect of the treatment applied to the real-world entity, and the method further comprises actioning the treatment on the real-world entity in dependence on the estimated treatment effect. Statement 14A. The method of Statement 13A, wherein one of: the real-world entity comprises a living being, and the treatment comprises a medical treatment to the living being or an environment thereof; or the real-world entity comprises a mechanical, electrical or electronic device or system, or an environment thereof; and the treatment comprises an act of maintaining, debugging, upgrading or controlling the device or system, or controlling the environment thereof; or the real-world entity comprises a network or software, and the treatment comprises an act of controlling the network or software. Statement 15A. The method of any of Statements 2A to 14A, wherein the inference network comprises a respective constituent inference network for each pair of variables in the feature vector, and e) comprises: for each identified confounded pair, inputting the input value of the respective parent variable and the input value of the selected variable into the respective inference network for the respective confounded pair, thereby resulting in the respective latent value of the respective confounder. Alternatively the inference network may comprise a common inference network operable to encode all the variables of the feature vector together into the respective latent for each pair. Statement 16A. The method of any of Statements 2A to 15A, wherein each of the first and second encoders and each of the decoders comprises a neural network. Statement 17A. The method of any of Statements 1A to 16A, wherein the inference network comprises a neural network. Statement 18A. The method of Statement 9A or any subsequent Statement when dependent thereon, wherein the method is performed on a server system of a first party, the server system comprising one or more server units at one or more sites; and the method further comprises, by the server system of the first party: providing an application programming interface, API, enabling a second party to contact the server system via a network; receiving a request from the second party over the network via the API; in response to the request, determining the estimated treatment effect on the target variable; and returning the estimated treatment effect to the second party over the network via the API. Statement 19A. A computer program embodied on computer-readable storage, wherein the computer program comprises a machine learning model comprising a plurality of first encoders, a plurality of second encoders, a decoder and an inference network; and wherein the computer- program further comprises instructions configured so as when run on one or more processors to perform the method of any of Statements 1A to 18A. Statement 20A. A system comprising: processing apparatus comprising one or more processors; and memory comprising one or more memory units, wherein the memory stores: a machine learning model comprising the first encoders, second encoders, decoders, and inference network of any preceding Statement; and code arranged to run on the processing apparatus and being configured so as when run to perform the method of any of Statements 1A to 18A. The first aspect (Statement 1) or any embodiment thereof may be used independently or in conjunction with the second aspect (Statement 1A) or any embodiment thereof. In some examples, the model 104/104’/104’’ and associated computer executable instructions are provided using any computer-readable media that are accessible by the computing equipment 102. Computer-readable media include, for example, computer storage media such as memory and communications media. Computer storage media include volatile and non-volatile, removable, and non-removable media implemented in any method or technology for storage of information such as computer readable instructions, data structures, program modules or the like. Computer storage media include, but are not limited to, Random Access Memory (RAM), Read-Only Memory (ROM), Erasable Programmable Read-Only Memory (EPROM), Electrically Erasable Programmable Read-Only Memory (EEPROM), persistent memory, phase change memory, flash memory or other memory technology, Compact Disk Read-Only Memory (CD-ROM), digital versatile disks (DVD) or other optical storage, magnetic cassettes, magnetic tape, magnetic disk storage, shingled disk storage or other magnetic storage devices, or any other non-transmission medium that can be used to store information for access by a computing apparatus.In contrast, communication media may embody computer readable instructions, data structures, program modules, or the like in a modulated data signal, such as a carrier wave, or other transport mechanism. As defined herein, computer storage media do not include communication media. Therefore, a computer storage medium should not be interpreted to be a propagating signal per se. Propagated signals per se are not examples of computer storage media. Although the computer storage medium may be described within the computing equipment 102, it will be appreciated by a person skilled in the art, that, in some examples, the storage is distributed or located remotely and accessed via a network or other communication link (e.g., using a communication interface). Other variants or use cases of the disclosed techniques may become apparent to the person skilled in the art once given the disclosure herein. The scope of the disclosure is not limited by the described embodiments but only by the accompanying claims.

Claims

CLAIMS 1. A computer-implemented method comprising: A) selecting a selected variable from among variables of a feature vector; B) sampling a temporal causal graph from a temporal graph distribution, the temporal graph distribution specifying probabilities of directed causal edges between different ones of the variables of the feature vector at a present time step, and from one of the variables of the feature vector at a preceding time step to one of the variables of the feature vector at the present time step; C) from among of the variables of the feature vector, identifying a present parent which is a cause of the selected variable in the present time step according to the temporal causal graph, and identifying a preceding parent which is a cause of the selected variable from the preceding time step according to the temporal causal graph; D) inputting an input value of each of the identified present and preceding parent into a respective encoder, resulting in a respective embedding of each of the present and preceding parents; E) combining the embeddings of the present and preceding parents, resulting in a combined embedding; and F) inputting the combined embedding into a decoder associated with the selected variable, resulting in a reconstructed value of the selected variable.

2. The method of claim 1, wherein: the temporal causal graph distribution sampled in B) specifies probabilities of directed causal edges existing from each of a plurality of variables of the feature vector from one or more preceding time steps, each to a variable of the feature vector at the present time step; and C) comprises, from among of the variables of the feature vector, identifying each present parent that is a cause of the selected variable in the present time step according to the temporal causal graph, and identifying each preceding parent variable which is a cause of the selected viable from a preceding time step according to the temporal causal graph.

3. The method of claim 2, further comprising: generating a history dependent noise term for the selected variable based on embeddings of the preceding parents; and combining the history dependent noise term with the reconstructed value of the selected variable, resulting in a simulated value of the reconstructed variable.

4. The method of claim 3, wherein embeddings of the present parents are not used to generate the history dependent noise term for the selected variable.

5. The method of any of claims 2 to 4, comprising: i) for a given training data point comprising a given combination of input values of the variables of the feature vector at the present times step and each preceding time step, repeating A)-F) over multiple selections, each selection selecting a different one of the variables of the feature vector as the selected variable thereby resulting in a respective reconstructed value, the multiple selections together thereby resulting in a reconstructed version of the training data point comprising the reconstructed values for the training data point; ii) evaluating a measure of difference between the training data point and the reconstructed version; and iii) training model parameters of the encoders, decoders, and temporal graph distributions, based on the evaluated measure; wherein the method comprises repeating i)-iii) over multiple input data points, each comprising a different combination of input values of the variables of the feature vector.

6. The method of claim 5, comprising repeating i)-iii) over multiple input data points, each comprising a different combination of input values of the variables of the feature vector.

7. The method of claim 5 or 6, wherein: each selection in i) further comprises a history dependent noise model generating a respective history dependent noise term based on embeddings of the preceding parents; wherein the training further comprises training model parameters the history dependent noise model based on the evaluated measure.

8. The method of any of claims 2 to 7, comprising: I) setting the input value of an intervened-on variable of the feature vector, other than the selected variable, to a specified value; and II) estimating an effect of the intervened-on variable on the selected variable based on the reconstructed value of the selected variable.

9. The method of claim 8, wherein: I) comprises setting the input value of a plurality of intervened-on variables of the feature vector, other than the selected variable, to respective specified values; and II) comprises estimating the effect of the plurality of intervened-on variables based on the reconstructed value of the selected variable.

10. The method of claim 9, wherein: I) comprises repeating A)-F) over multiple rounds, each using the same selected variable and setting the intervened-on variable to the same specified value, but performing the sampling of the temporal causal graph afresh each round, thereby resulting in a respective reconstructed value of the selected variable from each round; and II) comprises determining an average treatment effect of the intervened-on variable on the selected variable averaged over the multiple rounds.

11. The method of claim 10, wherein: each round sets a plurality of intervened-on variables and/or observed variables of the feature vector to specified values; and each round comprises an interior loop of D)-F) repeated around multiple iterations with the same sampled graph but for incrementing values of the present time step with each iteration, wherein each but the last iteration comprises, in addition to the target variable, reconstructing or simulating values of one or more further variables of the feature vector, and feeding back the simulated or reconstructed values as the specified input values of the next iteration.

12. The method of any of claims 8 to 11, wherein the intervened-on variable models a treatment on a real-world entity or an environment thereof and the target variable models an effect of the treatment applied to the real-world entity, and the method further comprises actioning the treatment on the real-world entity in dependence on the estimated treatment effect.

13. The method of claim 10 or 11, or claim 12 when dependent on claim 11 or 12, wherein each round further comprises: generating a history dependent noise term based on the embeddings of the preceding parents; and combining the history dependent noise term with the reconstructed value of the selected variable, resulting in a simulated value of the reconstructed variable; wherein said estimating based on the reconstructed values in II) comprises: estimating the effect of the intervened-on variable on the selected variable based on the simulated values.

14. A system comprising: processing apparatus comprising one or more processors; and memory comprising one or more memory units, wherein the memory stores: - a machine learning model comprising a plurality of encoders and a plurality of decoders, and - code arranged to run on the processing apparatus and being configured so as when run to perform a method comprising: A) selecting a selected variable from among variables of a feature vector; B) sampling a temporal causal graph from a temporal graph distribution, the temporal graph distribution specifying probabilities of directed causal edges between different ones of the variables of the feature vector at a present time step, and probabilities of causal edges existing from each of a plurality of variables of the feature vector at one or more preceding time steps to one or more variables of the feature vector at the present time step; C) from among of the variables of the feature vector, identifying each present parent which is a cause of the selected variable in the present time step according to the sampled temporal causal graph, and identifying each preceding parent which is a cause of the selected variable from any preceding time step according to the sampled temporal causal graph; D) inputting an input value of each of the identified present and preceding parents into a respective one of the encoders, resulting in a respective embedding of each of the present and preceding parents; E) combining the embeddings of the present and preceding parents, resulting in a combined embedding; and F) inputting the combined embedding into a respective one of the decoders associated with the selected variable, resulting in a reconstructed value of the selected variable.

15. A computer program embodied on non-transitory computer-readable storage, wherein the computer program comprises instructions configured so as when run on one or more processors to perform a method comprising: A) selecting a selected variable from among variables of a feature vector; B) sampling a temporal causal graph from a temporal graph distribution, the temporal graph distribution specifying probabilities of directed causal edges between different ones of the variables of the feature vector at a present time step, and probabilities of causal edges existing from each of a plurality of variables of the feature vector at one or more preceding time steps to one or more variables of the feature vector at the present time step; C) from among of the variables of the feature vector, identifying each present parent which is a cause of the selected variable in the present time step according to the sampled temporal causal graph, and identifying each preceding parent which is a cause of the selected variable from any preceding time step according to the sampled temporal causal graph; D) inputting an input value of each of the identified present and preceding parents into a respective first encoder, resulting in a respective embedding of each of the present and preceding parents; E) combining the embeddings of the present and preceding parents, resulting in a combined embedding; F) inputting the combined embedding into a decoder associated with the selected variable, resulting in a reconstructed value of the selected variable; and G) generating a history dependent noise term based on embeddings of the preceding parents, and combining the history dependent noise term with the reconstructed value of the selected variable, resulting in a simulated value of the reconstructed variable.