text:
Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

# Fast Maximum Likelihood Estimation via Equilibrium Expectation for Large Network Data

## Abstract

A major line of contemporary research on complex networks is based on the development of statistical models that specify the local motifs associated with macro-structural properties observed in actual networks. This statistical approach becomes increasingly problematic as network size increases. In the context of current research on efficient estimation of models for large network data sets, we propose a fast algorithm for maximum likelihood estimation (MLE) that affords a significant increase in the size of networks amenable to direct empirical analysis. The algorithm we propose in this paper relies on properties of Markov chains at equilibrium, and for this reason it is called equilibrium expectation (EE). We demonstrate the performance of the EE algorithm in the context of exponential random graph models (ERGMs) a family of statistical models commonly used in empirical research based on network data observed at a single period in time. Thus far, the lack of efficient computational strategies has limited the empirical scope of ERGMs to relatively small networks with a few thousand nodes. The approach we propose allows a dramatic increase in the size of networks that may be analyzed using ERGMs. This is illustrated in an analysis of several biological networks and one social network with 104,103 nodes.

## Introduction

Developing efficient approaches to data analysis and statistical inference is becoming increasingly important due to the widespread availability of large data sets in many fields of science. This is particularly the case for relational data typically taking the form of square arrays recording the presence of one or more relations among units of interest represented as network nodes1,2. Network representation of data is a fundamental tool for understanding and modeling a wide range of complex systems. Among the available models for relational data, exponential random graph models (ERGMs) are generally considered as “The most promising class of statistical models for expressing structural properties of social networks observed at one moment in time”3. ERGMs have found broad application in the analysis of social networks4, as well as biological networks5.

Formally, ERGMs may be written in the following form:

$$\pi (x,{\boldsymbol{\theta }})=\frac{1}{k({\boldsymbol{\theta }})}\exp \,(\sum _{A}{\theta }_{A}{z}_{A}(x)),$$
(1)

which expresses the probability of observing a network with a fixed number nodes in a given state x. Here zA(x) are network statistics from the space of states x, which are counts of theoretically meaningful or empirically relevant network subgraphs (often called “configurations” in ERGM terminology). The summation is over all configurations A, θA denotes the model parameter associated to zA(x), θ is the vector of all these parameters and k(θ) the normalizing constant ensuring that the probability distribution sums to one. The reader may also recognize equation (1) as an exponential family in canonical form (see Ch. 8 of Barndorff-Nielsen6), a Gibbs distribution, or a Markov random field7.

Different structural features are present in different networks. Yet, empirical network data are also characterized by recurrent structural regularities whose identification is crucial for understanding the behavior of complex network systems. For instance, much recent research has concerned “motifs,” small subgraphs occurring more frequently than might be expected by chance (and hence often similar to ERGMs configurations)8,9. Motifs have been considered as the building blocks of complex networks. To determine if a given motif is over-represented, the frequency of the motif in an observed network is typically compared to the average frequency in an appropriate random network (null model)10,11,12. ERGMs permit inference on the under or over-representation of specific configurations conditional on the presence of other configurations, or structural features. In ERGMs, network statistics may be defined for particular subgraphs that may be of general or contingent interest in the particular network under study like, for example, reciprocity or triadic closure3,4,13,14,15. Network nodes may have attributes16, the expression for zA(x) may be rather complicated, and the number of parameters may be large. Common examples of network statistics adopted in empirical research may be found elsewhere3,4,17.

The model parameters fit the observed network xobs if the following method of moments (MoM) condition18 is satisfied for all A:

$${E}_{\pi (\theta )}({z}_{A}(x))={z}_{A}({x}_{obs}),$$
(2)

where $${E}_{\pi (\theta )}({z}_{A}(x))={\sum }_{x}{z}_{A}(x)\pi (x,{\boldsymbol{\theta }})$$ denotes the expected network statistics with respect to the probability distribution (1) with parameters θ to be estimated; zA(xobs) denotes network statistics in the observed network. If the estimated θA is significantly larger than zero, then the corresponding network statistic zA occurs more frequently than might be expected by chance given all the other parameters of the model, and the corresponding configuration is over-represented. A well-known property of the exponential family (e.g., Ch. 8 of Barndorff-Nielsen6) is that Eπ(θ)(zA(x)) is a monotonically increasing function of θA. Thus, the estimated θA measures the corresponding structural features in observed networks.

The problem of parameter estimation (2) coincides, in our context, with the problem of maximum likelihood estimation (MLE). The computational challenge involved in using ERGMs is the intractable normalizing constant k(θ), that makes MLE computable only by Monte Carlo techniques14. Throughout the paper, we will be assuming that a MLE exists and that the model is non-degenerate14,19. It is known that this latter condition does not hold for all ERGM models (e.g., the edge-triangle model and other simple Markov random graph models13,19,20 with phase transitions, where non-convergent estimations are common), but more robust specifications, specifically the “alternating” statistics used in this paper (see the beginning of Section Results for the exact ERGM specifications adopted), have been shown to be much better behaved in terms of stability of statistics across a wide range of parameter values3.

Existing computational methods for MLE/MoM of ERGM parameters via equation (2) do not scale up easily to large data. Even though, to date, computational costs have constrained the scope of MLE, it remains widely adopted in numerous research settings, including the analysis of temporal networks21,22. Pseudo-likelihood and quasi-likelihood methods have been used when MLE cannot be computed, but it has been shown that these methods do not produce reliable results20,23,24. Estimation via conditional independence sampling methods, particularly snowball sampling, has been recently introduced to alleviate some of these issues of scale25,26,27.

Estimates of the model parameters θ may be achieved using a number of computational approaches such as, for example, Markov chain Monte Carlo maximum likelihood estimation (MCMCMLE)23,28,29,30,31, variants of stochastic approximation for the method of moments18,22,32, and Bayesian estimation33,34,35. Of these different approaches, Bayesian estimation of ERGMs is the most computationally intensive due to the double intractability of the posterior distribution34. In other types of models, however, new methods based on variational Bayesian inference36 can be relatively efficient, leveraging the sparse structure of large networks, e.g., graph partitioning (via blockmodeling analysis and/or analysis of modularity)37,38. These methods, however, have yet to be adopted for the estimation of ERGMs. In this paper, we derive an efficient Monte Carlo approach for MLE. MCMCMLE methods are somewhat faster than the MoM because they use importance sampling but are crucially dependent on the choice of the importance distribution. Stochastic approximation for the MoM is often more robust. For the network models we describe here, the existing MCMCMLE and MoM methods are similar in practical terms with respect to the network size that can be estimated. A considerable variety of stochastic approximation methods have been developed for different problems. In particular, efficient stochastic gradient methods39 are often adopted for maximum likelihood estimation when many independent observations are available. The MoM18 adopts the Robbins-Monro algorithm40 and Polyak-Ruppert averaging41 and is often used to solve (2) when the left side of (2) can be computed only by MCMC simulation42. In this paper, we compare the performance of our new approach with that of the MoM18.

MCMC simulation and, in particular, the Metropolis-Hastings algorithm43,44 may be used to generate a network x for fixed values of the parameters θ, which we denote by x(θ). Equation (2) formulates the inverse problem, i.e., to find the value of θ that, given xobs, satisfy (2), which we denote by θ(xobs). Computing θ(xobs) is much more computationally expensive than computing x(θ). The largest networks for which these methods have been applied to find the MLE of ERGM parameters contain at most a few thousand nodes. In contrast to existing approaches for MLE, what we propose does not require the simulation of a large number of Markov chains until convergence. Rather, it relies on properties of the Markov chain at equilibrium. For this reason, we call the proposed approach equilibrium expectation (EE).

The remainder of this paper is organized as follows: in the next section we give a brief description of Markov chain Monte Carlo, and then propose a new method, Equilibrium Expectation, for MLE of ERGM parameters. Then in the Results section, we demonstrate the performance of the EE algorithm on some simulated and empirical networks, of sizes well beyond what is currently possible in practice.

## Markov Chain Monte Carlo

If the value of the normalizing constant k(θ) in equation (1) is not known, the values of probability distribution (1) cannot be computed. Markov chain Monte Carlo simulation is typically used to address this issue. MCMC simulation allows approximation of the target probability (1) and computation of expected properties $${E}_{\pi (\theta )}({z}_{A}(x))={\sum }_{A}{z}_{A}(x)\pi (x,{\boldsymbol{\theta }})$$ of the model. The Metropolis-Hastings algorithm uses a Markov process that asymptotically reaches a unique stationary distribution π(θ). Given that the system is in state x the new state x′ is proposed with probability q(x → x′). A Markov process is uniquely defined by its transition probabilities, P(x → x′, θ) = q(x → x′)α(x → x′, θ), i.e., the probability of transitioning from any given state x to any state x′. Detailed balance with respect to a given distribution π(x, θ), which in turn implies that π(x, θ) is a stationary distribution for the Markov chain, is satisfied if the acceptance probability of the new state is given by

$$\alpha (x\to x^{\prime} ,{\boldsymbol{\theta }})=\,{\rm{\min }}\{\mathrm{1,}\,\frac{q(x^{\prime} \to x)\pi (x^{\prime} ,{\boldsymbol{\theta }})}{q(x\to x^{\prime} )\pi (x,{\boldsymbol{\theta }})}\}$$
(3)

Formally, the Metropolis-Hastings algorithm may be written as Algorithm 1.

If t is larger than the burn-in time, tB, of the Markov process, then xt values are drawn from π(x, θ), where tB is the time taken by the Markov chain to forget its initial state and reach the stationary regime. The algorithm generates zA(t) = zA(xt) sequences. From the Markov chain ergodic theorem it follows that, under regularity conditions, if the number of steps T is large, then the expected values of the model statistics Eπ(θ) (zA(x)) can be estimated by the average of zA(t) along the path of the Markov chain, $$\overline{{z}_{A}(t > {t}_{B})}$$44. The Metropolis-Hastings algorithm provides a general framework that results in a large number of different MCMC samplers that differ in the proposal q(x → x′). Verification of MCMC convergence is not a simple matter, and different methods have been suggested45. One heuristic rule for evaluating convergence of the Markov process xt, in our ERGM context, is to monitor the zA(xt) sequences.

Consider the expected change Eπ(θ)zA) of zA, Eπ(θ)zA) = Eπ(θ)(zA(xt+1) − zA(xt)). If equilibrium stationary distribution is reached, then statistics zA(xt) converge and fluctuate around their expected values Eπ(θ) (zA), which do not depend on t, and hence

$${E}_{\pi ({\boldsymbol{\theta }})}({\rm{\Delta }}{z}_{A})=0$$
(4)

We can write Eπ(θ)zA) as a function of the transition probabilities. Given xt = x, the expected change along the path of the Markov chain in zA(xt+1) − zA(xt) is obtained as

$${\rm{\Delta }}{z}_{A}(x,{\boldsymbol{\theta }})=\sum _{x^{\prime} }P(x\to x^{\prime} ,{\boldsymbol{\theta }})({z}_{A}(x^{\prime} )-{z}_{A}(x))$$
(5)

and the expected value of ΔzA(x, θ) with respect to π(x, θ) is

$${E}_{\pi ({\boldsymbol{\theta }})}({\rm{\Delta }}{z}_{A})=\sum _{x,x^{\prime} }\pi (x,{\boldsymbol{\theta }})P(x\to x^{\prime} ,{\boldsymbol{\theta }})({z}_{A}(x^{\prime} )-{z}_{A}(x))$$
(6)

From (4), it follows that after the MCMC burn-in, when equilibrium stationary distribution is reached, we have that for all the statistics zA(x),

$$\sum _{x,x^{\prime} }\pi (x,{\boldsymbol{\theta }})P(x\to x^{\prime} ,{\boldsymbol{\theta }})({z}_{A}(x^{\prime} )-{z}_{A}(x))=0$$
(7)

## Estimation strategy: Equilibrium Expectation

In the context of ERGMs, MLE of parameters is obtained from equation (2). Existing computational methods for MLE, such as MCMCMLE, the MoM, or Bayesian estimation, use iterative algorithms that successively modify θ until (2) is satisfied within given criteria for MCMC convergence. To this end, MCMC simulation is performed by drawing a large number of simulated networks for various values of θ, which we denote by xS(θ). The simulated network xS is a network drawn from probability distribution π(x, θ).

Each time a simulated network xS(θ) is drawn, the convergence criterion (7) should be satisfied, but this makes the standard estimation algorithms very computationally expensive. We suggest a much less computationally expensive approach for MLE.

We use equation (7) and rewrite it as:

$${E}_{\pi ({\boldsymbol{\theta }})}({\rm{\Delta }}{z}_{A}(x,{\boldsymbol{\theta }}))=0$$
(8)

Equation (8) reads as follows: if network x is drawn from probability distribution π(x, θ), then the expected value of ΔzA(x, θ) is zero. Equation (8) is valid only at equilibrium, that is, when π(x, θ) is the limiting equilibrium distribution of the Markov chain that has transition probabilities P(x → x′, θ). ΔzA(x, θ) may be computed via Monte Carlo integration42,46,47, as suggested in Section Estimation algorithm. The expectation with respect to π(x, θ) could be computed if a large Monte Carlo sample of networks independent and identically distributed (i.i.d.) from π(x, θ) were available i.e., $${x}_{{S}_{1}},{x}_{{S}_{2}},\cdots ,{x}_{{S}_{n}}$$. The existence of this large i.i.d. sample of networks is assumed here for expository purposes in deriving an estimation strategy, and later we will remove this assumption to develop an estimation method which can be applied to observed empirical networks. Making this assumption, the LHS of (8) may be computed by Monte Carlo integration as $${E}_{\pi ({\boldsymbol{\theta }})}({\rm{\Delta }}{z}_{A}(x{\boldsymbol{\theta }}))=\frac{1}{n}{\sum }_{i}{\rm{\Delta }}{z}_{A}({x}_{{S}_{i}},{\boldsymbol{\theta }})$$, and (8) may be approximated by

$$\frac{1}{n}\sum _{i}{\rm{\Delta }}{z}_{A}({x}_{{S}_{i}},{\boldsymbol{\theta }})=0$$
(9)

If we have a large sample of networks i.i.d. from π(x, θ), then we can efficiently compute the LHS of (9) and solve system of equations (9) with respect to θ. Thus, we can estimate θ. When MCMC simulation is performed, the number of steps should be larger than the burn-in time, which may be large. In contrast, ΔzA(x, θ) in (9) may be computed by Monte Carlo integration, there is no burn-in, and the number of steps may be small.

Equation (2) may be written as fA(xobs, θ) = 0, where fA(xobs, θ) = Eπ(θ) (zA(x) − zA(xobs)); hence, the true parameter values, $${{\boldsymbol{\theta }}}^{\ast }$$, may be estimated from

$$\frac{1}{n}\sum _{i}\,{f}_{A}({x}_{{S}_{i}},{{\boldsymbol{\theta }}}^{{\boldsymbol{\ast }}})=0$$
(10)

Typically, to estimate $${{\boldsymbol{\theta }}}^{\ast }$$, the estimation of each observed network is performed, and the resulting estimates $${\boldsymbol{\theta }}({x}_{{S}_{i}})$$ are averaged over the observations as follows: $$\hat{{\boldsymbol{\theta }}}$$* $$=\frac{1}{n}{\sum }_{i}{\boldsymbol{\theta }}({x}_{{S}_{i}})$$. If the network xS is very large, then $${\boldsymbol{\theta }}({x}_{{S}_{i}})$$ is the desired estimate of the true $${{\boldsymbol{\theta }}}^{\ast }$$. Thus, if the network xS is very large, then the true $${{\boldsymbol{\theta }}}^{\ast }$$ may be estimated from fA(xobs, θ) = 0, and the summation in equation (10) may be dropped. In statistics and statistical physics, this property is called the ergodicity of systems: the time averaging is equivalent to the ensemble averaging. Here, time averaging is the average over networks generated by the Markov process, and ensemble averaging is the average over the space of all the system’s states, which grows with the network size. If the system is ergodic, then this property may also be applied to equation (9) i.e., if the network xs is very large, then the summation in equation (9) may be dropped. Thus, the true $${{\boldsymbol{\theta }}}^{\ast }$$ value may be found from the following equilibrium expectation condition: for all A,

$${\rm{\Delta }}{z}_{A}({x}_{S},{{\boldsymbol{\theta }}}^{{EE}})=\mathrm{0,}$$
(11)

where ΔzA(xS, θ) is given by (5). We can find the value of θEE that, given xS, satisfies (11); this value is denoted as θEE(xS). If the network xS is very large, then θEE(xS) is the desired estimate of $${{\boldsymbol{\theta }}}^{\ast }$$. Otherwise, a large sample of networks is required as described above, and the true $${{\boldsymbol{\theta }}}^{\ast }$$ value may be found from (9).

The estimation method described above may be theoretically interesting, but is of limited value for the estimation of real-world data, for which we have only a single observed network xobs drawn from some unknown probability distribution, rather than a large i.i.d. Monte Carlo sample of networks drawn from a known probability distribution, as we assumed earlier.

It is not, however, lacking in usefulness when applied to such an observed network. In fact, it is actually Contrastive Divergence48 (CD) as applied to ERGM parameter estimation31,49,50,51. The CD algorithm, rather than running an MCMC simulation to convergence, instead starts from the observed data and makes only some number k of MCMC updates (the CD algorithm with k updates is called CD-k). CD-1 has been shown to be equivalent to maximum pseudo-likelihood under certain conditions52 and it has also been shown that CD-k forms a series of increasingly close approximations to the MLE as k → ∞49.

To see that the estimation strategy we have just described is equivalent (but independently developed) to CD-1, it may be useful to consider the implementation described by Krivitsky31 of using CD to find the initial estimates to seed MCMCMLE. In this implementation, an ERGM MCMCMLE implementation is simply converted to CD-k by modifying the MCMC sampler to start from xobs and reverting the chain to xobs every k steps. And so if k = 1 then the network is not actually modified and the CD-1 estimate of θ is just the solution of (10). Details of the CD-1 algorithm as applied in the work described here are given in the Supplementary Information.

As previously noted, maximum pseudo-likelihood (and hence CD-1) does not necessarily produce reliable results for ERGM estimation. However they (and CD-k for k > 1) are useful for finding initial θ values to seed methods such as MCMCMLE that are sensitive to initialization conditions31,53. In this paper we indeed use CD-1 to find initial parameter estimates.

We will now describe how this strategy can be modified, for more realistic applications, having an observed network drawn from some unknown probability distribution. Let $${{\boldsymbol{\pi }}}^{\ast }$$(x) be the unknown probability distribution from which the real-world data xobs are drawn. π(x, θ) is then the probability distribution corresponding to the model specification we have chosen (in this paper the ERGM model for some specified choices of network statistics). We assume this model to be appropriate for the observed data – we do not discuss model selection here – and, as discussed in the Introduction, we also assume the MLE exists and the model is not degenerate. Formally, we can write that for any θ, $${{\pi }}^{\ast }$$(x) ≠ π(x, θ). Equation (7) can be applied only if the equilibrium stationary distribution of the Markov chain follows a statistical model π(x, θ). If $${{\pi }}^{\ast }$$(x) ≠ π(x, θ), then the conditions of Equation (7) are not satisfied, and it cannot be applied. In other words, if observed networks are not drawn from π(x, θ), then the LHS of equation (7) is not equal to the LHS of equation (9). In this case, we cannot compute the LHS of (7), and hence we cannot find a value of θ such that (7) is satisfied. However, if a solution to (2) exists, then it is possible to find xS drawn from π(x, θ) such that the following condition is satisfied:

$${z}_{A}({x}_{S})={z}_{A}({x}_{obs})\mathrm{.}$$
(12)

For network xS, drawn from probability distribution π(x, θ), (11) may be used as described above. Using (10) and (12), the solution of (2) may be obtained for an observed network drawn from an unknown probability distribution. The network drawn from π(x, θ) may be obtained by an MCMC simulation. In the estimation strategy just described, the parameters θ are iteratively adjusted until (11) is satisfied. To estimate parameters for networks drawn from an unknown distribution $${{\pi }}^{\ast }$$(x), both θ and x are modified: θ is modified in the same way as before, and x is also modified (equilibrated) so that, after some burn-in, x will be drawn from the probability distribution π(x, θ).

Indeed, starting from the previously described method to find the solution of (11), we can derive the EE algorithm (described in the following section) for the solution of (11), (12). The EE algorithm for the estimation of an observed network xobs performs Metropolis-Hastings moves, that is, it actually makes changes to a network for accepted proposals, a step which was not necessary for solving (11). The EE algorithm requires MCMC simulation, but – in contrast to MCMCMLE, the MoM or Bayesian estimation – the EE algorithm does not draw many simulated networks for various θ values and, therefore, is considerably faster.

## Estimation algorithm

The EE algorithm is described in Algorithm 2. To provide some intuition, the description is maintained at a rather general level here. A more detailed and formal description of the EE algorithm is included in the Supplementary Information. Let m = 1000 and M indicate the number of steps of the EE algorithm. KA is a positive constant, the values of which could be for example KA ≈ 10−4 (∂ΔzA(xobs, θ)/∂θA)−2; a better choice is suggested in the Supplementary Information. The final parameter estimate from the CD-1 algorithm is used as the starting point for the EE algorithm, but the efficiency of the EE algorithm is more sensitive to KA values than to the initial parameter estimates θ0.

Steps 4–5 of the EE algorithm are identical to steps 3–4 of the Metropolis-Hastings algorithm presented in Section Chain Monte Carlo: the move x → x′ is proposed, and the acceptance probability (3) is calculated. If the move to x′ is accepted, then the change in the sufficient statistics zA(x′) − zA(x) is computed.

Step 13 modifies θ iteratively until (11) is satisfied. This parameter update depends on the property of the exponential family that Eπ(θ)(zA(x)) is a monotonically increasing function of θA, which was mentioned in the Introduction. Because of this, for each configuration A we want to increase θA when dzA is negative and decrease it when dzA is positive, in order to find θA such that dzA is zero. This is achieved by adding $$-{K}_{A}\cdot {\rm{sign}}(d{z}_{A})d{z}_{A}^{2}$$ to the current θA value (recall that KA is positive).

Different root-finding algorithms may be applied to solve the system of equation (11) when ΔzA(x, θ) is available. For example, bisection or quasi-Newton methods46,54 may be applied if ΔzA(x, θ) were a continuous functions of θ. Efficient root finding methods use derivative information. The precise derivative information is either not available or very computationally expensive, but it may be shown (see Supplementary Information) that, for any x,

$$\partial {\rm{\Delta }}{z}_{A}(x,{\boldsymbol{\theta }})/\partial {\theta }_{A}\ge 0.$$
(13)

It is trivial to find roots of continuous monotonic functions, and step 13 implements one possibility. A more robust approach to solving ΔzA(x, θ) = 0 would be to use stochastic approximation methods, which are typically applied to solve (2).

Since x is not constant, zA(x) is also not constant. Recently, an auxiliary parameter MCMC method was proposed to perform MCMC simulation, constraining the value of one of the network statistics. Only the network statistic zL(x), which is the edge count, was constrained to the observed value, and this was achieved by using a special proposal q(x → x′): the IFD sampler55. The corresponding parameter θL was adapted55 in the same way as done in step 13 of the EE algorithm. In the present paper, all the network statistics zA(x) are constrained and, crucially, without the need to develop a special proposal distribution. Indeed, they may be constrained by exploiting (13), i.e., the monotonic dependence of ΔzA(x, θ) and Eπ(θ)(zA(x)) on θA, and by properly modifying the θA values. We need an algorithm that converges to the values of θ and x that satisfy (11) and (12) for all the network statistics zA(x). This is achieved by accumulating the accepted change statistics (line 7). As a result, dzA = dzA(t) + ΔzA(x, θ) and dzA = zA(x) − zA(xobs). If equilibrium is reached, then x = xS is drawn from π(x, θ) and hence ΔzA(x, θ) = 0 and (11) is satisfied. In addition, if dzA = 0, then zA(x) = zA(xobs) and (12) is satisfied. It is shown in Section Estimation strategy: Equilibrium Expectation that the solution of (11) and (12) gives the solution of (2).

The EE algorithm produces θA(t) sequences. The number of steps M should be large enough for θA(t) to converge. After θA(t) reaches convergence, its values fluctuate around some constant values θA, which we estimate as $${\theta }_{A}=\overline{{\theta }_{A}(t > {t}_{B})}$$. In addition to θA(t), the sequences dzA(t) = zA(t) − zA(xobs) are generated by the EE algorithm. It is convenient to plot both zA − zA(xobs) and θA(t) sequences to visually check for convergence of the algorithm. We use the t-ratio test to check if (12) is satisfied:

$${\tau }_{A}=\frac{\overline{{z}_{A}(t > {t}_{B})}-{z}_{A}({x}_{obs})}{SD({z}_{A}(t > {t}_{B}))},$$
(14)

where tB is the burn-in time, starting from which convergence is observed from the plots of θA(t) and zA(t) − zA(xobs), and $$\overline{{z}_{A}(t > {t}_{B})}$$ denotes the network statistic averaged over t > tB, while SD() is the standard deviation operator. We propose an appropriate indication of good convergence so that the suggested algorithm gives a converged solution to (2) if, for all A, (i) |τA| < 0.1 and (ii) θA converges.

Note that if KA = 0, then the EE algorithm does not differ from the Metropolis-Hastings algorithm. MCMC estimation algorithms make use of MCMC simulation; hence, the efficiency of MCMC estimation algorithms depends on the efficiency of MCMC samplers. In turn, the efficiency of the Metropolis-Hastings algorithm depends on the proposal q(x → x′). The same holds also for the EE algorithm: the algorithm may be used with different proposals q(x → x′), and its efficiency depends on a good choice of the proposal distribution. Proposals for ERGMs are reviewed elsewhere55. Both the updating step of the MoM18 and that of the EE algorithm depend on zA(x) − zA(xobs). However, in the case of the MoM, x is the equilibrium network configuration, which can be obtained by, e.g., the Metropolis-Hastings algorithm if the MCMC simulation time is larger than the burn-in time of the Markov process. In the EE algorithm, x is a current non-equilibrium state. The MoM and MCMCMLE28 require many converged outputs of the Metropolis-Hastings algorithm, while the EE algorithm does not need such outputs. Instead, the EE algorithm generates one converged output.

## Results

First, we test the EE algorithm by computing the bias and the standard deviation of the estimates that this algorithm produces. Simulated networks of various sizes N (here, N is the number of nodes in a network) are generated by the Metropolis-Hastings algorithm (as implemented in the PNet56 program). The ERGM model is defined by the AS, AT, A2P, ρ and ρB network statistics. The network statistics we use in this paper were suggested by Snijders et al.3,20 and are detailed in the Supplementary Information. The true values of the corresponding model parameters $${\theta }_{A}^{\ast }$$ are represented by the horizontal lines in Fig. 1, with the estimated values of the parameters shown as boxplots, for both the EE algorithm and a variant of the stochastic approximation algorithm (MoM)18 for MLE, widely used for ERGM estimation. It is clear that both EE and the MoM give accurate estimates of the true $${\theta }_{A}^{\ast }$$ values. Both the variance and the bias $$|{\theta }_{A}-{\theta }_{A}^{\ast }|$$ of the estimates obtained using different methods have similar values. In Fig. 2, we report the computational times needed to obtain the corresponding estimates. As described in the Methods section, steps 4–6 of the EE algorithm are equivalent to steps 3–5 of the Metropolis-Hastings algorithm, used in the MoM. In our implementation of the MoM and EE algorithms, these steps are carried out in exactly the same way. This allows the evaluation of the efficiency of the MoM and EE algorithms by comparing the computational times for the corresponding estimations. Both the MoM and EE algorithms are used with the Basic sampler55, and estimation is performed on a Cray XC50 machine available at the Swiss National Supercomputing Centre (CSCS). Figure 2 shows that for small networks (N = 1,000), the efficiencies of the EE and MoM algorithms are not very different, while for larger networks (N = 10,000), the EE algorithm is 2 orders of magnitude faster than the MoM. These results suggest that the EE algorithm may be used for the analysis of large datasets.

For a more complete comparison of the MoM and EE, we now compare their respective performance on empirical datasets. Using the MoM, it is possible to estimate simulated networks with 10,000 nodes within several hours, but the estimation of empirical networks takes longer to converge. Despite several attempts to improve the efficiency of the computational methods for MLE, the existing algorithms are still too computationally expensive for empirical networks of this or larger size. Using the more efficient IFD sampler55, we can estimate six biological networks via both EE and the MoM. The estimation results are reported in the Supplementary Information. The results reported in Supplementary Tables S2 and S3 show that the MoM and EE produce equivalent estimates, but EE may be two orders of magnitude faster than the MoM, even for relatively small empirical networks with 1,781 and 5,038 nodes. Although we cannot estimate larger networks with the MoM, the results reported in Fig. 2 suggest that for larger networks, EE may be orders of magnitude faster.

For an illustrative example of a biological network, we study an Arabidopsis thaliana protein-protein interaction (PPI) network57,58, including some protein attributes on the nodes. In this network, nodes represent proteins, and (undirected) edges represent literature-curated binary interactions between proteins. Arabidopsis thaliana proteins are annotated with various properties, as described in the Supporting Online Material of Arabidopsis Interactome Mapping Consortium58. We make use of the following protein attributes. We consider one binary “plant-specific” node attribute (genes defined as plant-specific, absent from other eukaryotic lineages58) and define activity ρ and interaction ρB network statistics on this attribute. We consider 2 categorical attributes “E class” and “kinase-phosphorylated”, derived from the Supporting Online Material of Arabidopsis Interactome Mapping Consortium58 as described in the Supplementary Information, and define mismatch network statistics on these attributes. The estimated parameters of the resulting model are given in Table 1. The significant positive AT parameter indicates that triangle motifs are significantly over-represented. The positive plant-specific interaction parameter shows that plant-specific proteins interact preferentially with each other (with the negative activity effect showing that, relative to other proteins, they are less likely to interact at all). The positive mismatch kinase-phosphorylated effect shows the propensity of kinases to interact with phosphorylated proteins, as we would expect given the action of kinases to phosphorylate proteins. The positive mismatch E class shows the propensity of proteins to interact if they belong to different E classes59.

We also demonstrate the performance of the algorithm in a case study involving the analysis of a larger empirical social network. More specifically, we study the Livemocha network of friendship relations among members of an online language learning community available from the repository at http://konect.uni-koblenz.de60,61. Founded in 2007, Livemocha was an online community of peer-to-peer language learning that provided language instruction to approximately 12 million users from more than 190 countries. Livemocha closed in 2016 after being included by the Times magazine in its list of 50 best websites in 2010. The Livemocha dataset that we analyze consists of an undirected network with 104,103 nodes and 2,193,083 ties. We estimate AS, AT and edge parameters using the EE algorithm and IFD sampler. The convergence criteria |τA| < 0.1 is satisfied, and the values of estimated parameters θA(t) are presented in Fig. 3 as a function of the algorithm step t. One can see that θA(t) converge. The θA(t = 0) values are obtained using the CD-1 algorithm, and the smaller KA values of the EE algorithm are used for t > 2.5 × 106 (the computational details are given in the Supplementary Information). The values of the estimated parameters are significantly different from zero. The θA(t) values are averaged over the last 5 × 105 steps, and the following estimates of the model parameters are obtained: θL = −19.321, θAS = 2.7355, and θAT = 0.6453. A different convergence test was suggested by Snijders for the MoM algorithm: estimated parameter values may be used to estimate Eπ(θ)(zA(x)) in (2) via the Metropolis-Hastings algorithm, and the t-statistics18 may be computed. We checked that this convergence test is also satisfied: for all A, the absolute values of the t-statistics are less than 0.3. However, the convergence test, as suggested by Snijders, is computationally expensive. When the EE algorithm is applied, it is more convenient to check the values of τA and to visualize θA(t) dependencies.

## Discussion

We propose a Monte Carlo Markov chain based approach to the maximum likelihood estimation of parameters of statistical models with intractable normalizing constants belonging to the linear exponential family. The algorithm we propose is similar to the Metropolis-Hastings algorithm, but allows Monte Carlo simulation to be performed while constraining the values of model statistics so that the equilibrium expectation condition (11) is satisfied. In contrast to the widely adopted Metropolis-Hastings algorithm, which is guaranteed to converge within the limit of infinite simulation time, the EE algorithm provides no such guarantee. Further work is required to establish the general convergence conditions of the algorithm. This work is currently in progress.

We demonstrate the merits of the algorithm by estimating the parameters of ERGMs on large network data sets. ERGMs are popular statistical models for the analysis of single-observation network data. ERGMs may be used to determine motif significance, and are widely used to connect motifs to structural features in social and other kinds of networks. We compute the bias of the estimates produced by the EE algorithm and show that its value is close to the bias of a commonly adopted computational method for MLE. However, EE is faster, and the estimation time increases with the number of network nodes as N1.5. These results suggest that the EE algorithm may be relied upon to support statistical inference on very large complex systems with network-like local dependencies. We then study several biological networks and a large social network. We show that the triangle motif is significantly over-represented in all the networks studied. We show that accurate maximum likelihood estimates of ERGM parameters may be obtained for a large empirical network with 104,103 nodes and 2, 193,083 ties. The smallest network for which we applied the EE algorithm is a network with 418 nodes (see Supplementary Information).

The algorithm is inappropriate for the curved exponential family30 when the number of parameters to be estimated differs from the number of statistics. Currently, we are applying the EE algorithm to study directed networks for which different network statistics need to be computed3. Future research directions involve an attempt to couple the EE algorithm with the expectation-maximization algorithm62 for incomplete data, and an extension to bipartite networks. In addition, we believe the computational approach we have proposed in this paper may be applicable to a wider range of statistical models for data characterized by complex network-like dependencies.

### Data availability

The datasets analyzed in the current study are available in the following repositories: Livemocha http://konect.uni-koblenz.de/networks/livemocha, http://socialcomputing.asu.edu/datasets/Livemocha; A. thaliana PPI: http://interactome.dfci.harvard.edu/A_thaliana; Yeast PPI: igraph63 Nexus repository (as of March 2017, the igraph Nexus repository is no longer available; we downloaded this data set on 10 November 2016, and it is available upon request); Human PPI: https://icon.colorado.edu/#!/networks, http://interactome.dfci.harvard.edu/H_sapiens; C. elegans PPI: (downloaded from the web address specified in Huang et al.64 on 2 March 2017, which is no longer available, data available on request); E. coli regulatory: http://www.statnet.org/index.shtml; Drosophila optic medulla: https://icon.colorado.edu/#!/networks, http://openconnecto.me/graph-services/download/.

### Code availability

The EE algorithm is implemented in the Estimnet program, available from the authors upon request and from the webpage http://www.estimnet.org. The algorithm may be accessed from the webpage https://github.com/Byshkin/EquilibriumExpectation.

## References

1. Borgatti, S. P., Mehra, A., Brass, D. J. & Labianca, G. Network analysis in the social sciences. Science 323, 892–895 (2009).

2. Butts, C. T. Revisiting the foundations of network analysis. Science 325, 414–416 (2009).

3. Snijders, T. A. B., Pattison, P. E., Robins, G. L. & Handcock, M. S. New specifications for exponential random graph models. Sociol. Methodol. 36, 99–153 (2006).

4. Lusher, D., Koskinen, J. & Robins, G. Exponential random graph models for social networks: Theory, methods, and applications (Cambridge University Press, 2013).

5. Saul, Z. M. & Filkov, V. Exploring biological network structure using exponential random graph models. Bioinformatics 23, 2604–2611 (2007).

6. Barndorff-Nielsen, O. Information and exponential families in statistical theory (John Wiley & Sons, 2014).

7. Geman, S. & Geman, D. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. In Readings in Computer Vision, 564–584 (Elsevier, 1987).

8. Milo, R. et al. Network motifs: simple building blocks of complex networks. Science 298, 824–827 (2002).

9. Shen-Orr, S. S., Milo, R., Mangan, S. & Alon, U. Network motifs in the transcriptional regulation network of Escherichia coli. Nat. Genet. 31, 64–68 (2002).

10. Artzy-Randrup, Y., Fleishman, S. J., Ben-Tal, N. & Stone, L. Comment on “Network motifs: simple building blocks of complex networks” and “Superfamilies of evolved and designed networks”. Science 305, 1107 (2004).

11. Ciriello, G. & Guerra, C. A review on models and algorithms for motif discovery in protein–protein interaction networks. Brief. Funct. Genomic. Proteomic. 7, 147–156 (2008).

12. Kovanen, L., Kaski, K., Kertész, J. & Saramäki, J. Temporal motifs reveal homophily, gender-specific patterns, and group talk in call sequences. Proc. Natl. Acad. Sci. USA 110, 18070–18075 (2013).

13. Frank, O. & Strauss, D. Markov graphs. J. Am. Stat. Assoc. 81, 832–842 (1986).

14. Hunter, D. R., Krivitsky, P. N. & Schweinberger, M. Computational statistical methods for social network models. J. Comput. Graph. Stat. 21, 856–882 (2012).

15. Newman, M. E., Watts, D. J. & Strogatz, S. H. Random graph models of social networks. Proc. Natl. Acad. Sci. USA 99, 2566–2572 (2002).

16. Newman, M. E. & Clauset, A. Structure and inference in annotated networks. Nat. Commun. 7, 11863 (2016).

17. Pallotti, F., Lomi, A. & Mascia, D. From network ties to network structures: Exponential random graph models of interorganizational relations. Qual. Quant. 47, 1665–1685 (2013).

18. Snijders, T. A. B. Markov chain Monte Carlo estimation of exponential random graph models. J. Soc. Struct. 3, 1–40 (2002).

19. Handcock, M. S. Statistical models for social networks: Inference and degeneracy. In Dynamic Social Network Modeling and Analysis: Workshop Summary and Papers, 229–240 (National Academies Press, 2003).

20. Robins, G., Snijders, T. A. B., Wang, P., Handcock, M. & Pattison, P. Recent developments in exponential random graph (p*) models for social networks. Soc. Networks 29, 192–215 (2007).

21. Snijders, T. A. B., Koskinen, J. & Schweinberger, M. Maximum likelihood estimation for social network dynamics. Ann. Appl. Stat. 4, 567–588 (2010).

22. Snijders, T. A. B. The statistical evaluation of social network dynamics. Sociol. Methodol. 31, 361–395 (2001).

23. Hummel, R. M., Hunter, D. R. & Handcock, M. S. Improving simulation-based algorithms for fitting ERGMs. J. Comput. Graph. Stat. 21, 920–939 (2012).

24. van Duijn, M. A., Gile, K. J. & Handcock, M. S. A framework for the comparison of maximum pseudo-likelihood and maximum likelihood estimation of exponential family random graph models. Soc. Networks 31, 52–62 (2009).

25. Pattison, P. E., Robins, G. L., Snijders, T. A. B. & Wang, P. Conditional estimation of exponential random graph models from snowball sampling designs. J. Math. Psychol. 57, 284–296 (2013).

26. Stivala, A. D., Koskinen, J. H., Rolls, D. A., Wang, P. & Robins, G. L. Snowball sampling for estimating exponential random graph models for large networks. Soc. Networks 47, 167–188 (2016).

27. Thiemichen, S. & Kauermann, G. Stable exponential random graph models with non-parametric components for large dense networks. Soc. Networks 49, 67–80 (2017).

28. Geyer, C. J. & Thompson, E. A. Constrained Monte Carlo maximum likelihood for dependent data. J. Roy. Stat. Soc. B Met. 54, 657–699 (1992).

29. Geyer, C. J. Markov chain Monte Carlo maximum likelihood. In Keramides, E. M. (ed.) Computing Science and Statistics: Proceedings of the 23rd Symposium on the Interface, 156–163 (Interface Foundation of North America, 1991).

30. Hunter, D. R. & Handcock, M. S. Inference in curved exponential family models for networks. J. Comput. Graph. Stat. 15, 565–583 (2006).

31. Krivitsky, P. N. Using contrastive divergence to seed Monte Carlo MLE for exponential-family random graph models. Comput. Stat. Data Anal. 107, 149–161 (2017).

32. Okabayashi, S. et al. Long range search for maximum likelihood in exponential families. Electron. J. Stat. 6, 123–147 (2012).

33. Lehmann, E. L. & Casella, G. Theory of point estimation (Springer Science & Business Media, 2006).

34. Caimo, A. & Friel, N. Bayesian inference for exponential random graph models. Soc. Networks 33, 41–55 (2011).

35. Liang, F., Jin, I. H., Song, Q. & Liu, J. S. An adaptive exchange algorithm for sampling from distributions with intractable normalizing constants. J. Am. Stat. Assoc. 111, 377–393 (2016).

36. Blei, D. M., Kucukelbir, A. & McAuliffe, J. D. Variational inference: A review for statisticians. J. Am. Stat. Assoc. 112, 859–877 (2017).

37. Zhang, P. & Moore, C. Scalable detection of statistically significant communities and hierarchies, using message passing for modularity. Proc. Natl. Acad. Sci. USA 111, 18144–18149 (2014).

38. Zhang, P., Krzakala, F., Reichardt, J. & Zdeborová, L. Comparative study for inference of hidden classes in stochastic block models. J. Stat. Mech. Theory Exp. 2012, P12021 (2012).

39. Roux, N. L., Schmidt, M. & Bach, F. R. A stochastic gradient method with an exponential convergence rate for finite training sets. Adv. Neural Inf. Process. Syst. 25, 2663–2671 (2012).

40. Robbins, H. & Monro, S. A stochastic approximation method. Ann. Math. Stat. 40, 400–407 (1951).

41. Polyak, B. T. & Juditsky, A. B. Acceleration of stochastic approximation by averaging. SIAM J. Contr. Optim. 30, 838–855 (1992).

42. Christian, P. R. & Casella, G. Monte Carlo statistical methods (Springer, 1999).

43. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. & Teller, E. Equation of state calculations by fast computing machines. J. Chem. Phys. 21, 1087–1092 (1953).

44. Hastings, W. K. Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57, 97–109 (1970).

45. Cowles, M. K. & Carlin, B. P. Markov chain Monte Carlo convergence diagnostics: a comparative review. J. Am. Stat. Assoc. 91, 883–904 (1996).

46. Miranda, M. J. & Fackler, P. L. Applied computational economics and finance (MIT press, 2004).

47. Newman, M. & Barkema, G. Monte Carlo Methods in Statistical Physics chapter 1–4 (Oxford University Press: New York, USA, 1999).

48. Hinton, G. E. Training products of experts by minimizing contrastive divergence. Neural Comput. 14, 1771–1800 (2002).

49. Asuncion, A., Liu, Q., Ihler, A. & Smyth, P. Learning with blocks: Composite likelihood and contrastive divergence. In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, 33–40 (2010).

50. Hummel, R. M. Improving estimation for exponential-family random graph models. Ph.D. thesis, The Pennsylvania State University, https://etda.libraries.psu.edu/catalog/11493 (2010).

51. Fellows, I. E. Why (and when and how) contrastive divergence works. arXiv preprint arXiv:1405.0602 (2014).

52. Hyvärinen, A. Consistency of pseudolikelihood estimation of fully visible Boltzmann machines. Neural Comput. 18, 2283–2292 (2006).

53. Carreira-Perpiñan, M. A. & Hinton, G. E. On contrastive divergence learning. In Proceedings of the Tenth International Workshop on Artificial Intelligence and Statistics, 33–40 (2005).

54. Ortega, J. M. & Rheinboldt, W. C. Iterative solution of nonlinear equations in several variables (SIAM, 1970).

55. Byshkin, M. et al. Auxiliary parameter MCMC for exponential random graph models. J. Stat. Phys. 165, 740–754 (2016).

56. Wang, P., Robins, G. & Pattison, P. PNet: program for the estimation and simulation of p* exponential random graph models. Department of Psychology, The University of Melbourne (2009).

57. Swarbreck, D. et al. The Arabidopsis Information Resource (TAIR): gene structure and function annotation. Nucleic Acids Res. 36, D1009–D1014 (2008).

58. Arabidopsis Interactome Mapping Consortium. Evidence for network evolution in an Arabidopsis interactome map. Science 333, 601–607 (2011).

59. Mazzucotelli, E. et al. The E3 ubiquitin ligase gene family in plants: regulation by degradation. Curr. Genom. 7, 509–522 (2006).

60. Liaw, M.-L. Review of livemocha. Lang. Learn. Technol. 15, 36–40, http://www.lltjournal.org/item/2722 (2011).

61. Zafarani, R. & Liu, H. Social computing data repository at ASU, http://socialcomputing.asu.edu (2009).

62. Dempster, A. P., Laird, N. M. & Rubin, D. B. Maximum likelihood from incomplete data via the EM algorithm. J. Roy. Stat. Soc. B Met. 1–38 (1977).

63. Csárdi, G. & Nepusz, T. The igraph software package for complex network research. InterJournal Complex Systems, 1695, http://igraph.org (2006).

64. Huang, X.-T., Zhu, Y., Chan, L. L. H., Zhao, Z. & Yan, H. An integrative C. elegans protein–protein interaction network with reliability assessment based on a probabilistic graphical model. Mol. BioSyst. 12, 85–92 (2016).

65. Wickham, H. ggplot2: Elegant Graphics for Data Analysis, http://ggplot2.org (Springer-Verlag New York, 2009).

## Acknowledgements

The authors are grateful to Dr. Peng Wang for making the source code for the PNet program available. Mark Girolami, Illia Horenko, Ritabrata Dutta, Pavel Krivitsky and Akihiko Nishimura provided helpful discussions on earlier versions of the manuscript. We are grateful to Federica Bianchi for her careful and valuable assistance and advice on earlier drafts. We thank Melbourne Bioinformatics at the University of Melbourne, grant number VR0261, and the Swiss National Supercomputing Centre for high performance computing facilities. M.B. thanks Svitlana Byshkina for her patience. This work was supported by Swiss National Science Foundation grant number 167362, National Research Program 75 (“Big Data”) awarded to A. Lomi, and by Swiss National Science Foundation grant number 105218_163196 (“Statistical Inference on Large-Scale Mechanistic Network Models”) awarded to A. Mira.

## Author information

Authors

### Contributions

M.B., A.S., A.M., G.R. and A.L. conceived the research, designed the research, conducted the analyses and wrote the manuscript. M.B. conceived the original EE algorithm.

### Corresponding author

Correspondence to Alessandro Lomi.

## Ethics declarations

### Competing Interests

The authors declare no competing interests.

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Rights and permissions

Reprints and Permissions

Byshkin, M., Stivala, A., Mira, A. et al. Fast Maximum Likelihood Estimation via Equilibrium Expectation for Large Network Data. Sci Rep 8, 11509 (2018). https://doi.org/10.1038/s41598-018-29725-8

• Accepted:

• Published:

• DOI: https://doi.org/10.1038/s41598-018-29725-8

• ### Testing biological network motif significance with exponential random graph models

• Alex Stivala
• Alessandro Lomi

Applied Network Science (2021)

• ### Tie-formation process within the communities of the Japanese production network: application of an exponential random graph model

• Hazem Krichene
• Abhijit Chakraborty
• Masaaki Terai

Applied Network Science (2019)


text:
ISO TC 211 Geographic information/Geomatics

# ISO/TC 211 Geolexica

## Concept “coordinate conversion”

Term ID

2350

source

ISO 19111:2019, (E), 3.1.6

eng

### coordinate conversion

coordinate operation that changes coordinates in a source coordinate reference system to coordinates in a target coordinate reference system in which both coordinate reference systems are based on the same datum

Note to entry: A coordinate conversion uses parameters which have specified values.

EXAMPLE 1 A mapping of ellipsoidal coordinates to Cartesian coordinates using a map projection.

EXAMPLE 2 Change of units such as from radians to degrees or from feet to metres.

ORIGIN: ISO/TC 211 Glossary of Terms - English (last updated: 2020-06-02)

JSON

/api/concepts/2350.json

SKOS in JSON-LD

/api/concepts/2350.jsonld

SKOS in RDF

/api/concepts/2350.ttl

info

• status: valid
• classification: preferred
• date accepted: 2019-01-31

Review

last review performed:
(2019-01-31)
status:
final
decision:
accepted
decision event:
Publication of document ISO 19111:2019(E)
notes:
This entry supersedes the entry for coordinate conversion in ISO 19111:2007(E)

format: plain


text:
# How to specify the URL of the page the link goes to in HTML?

HTMLWeb DevelopmentFront End Technology

Use the href attribute to specify the URL of the page the link goes to in HTML.

## Example

You can try to run the following code to implement href attribute −

<!DOCTYPE html>
<html>
</html>


text:
What is the TF tree? Set frame_id to base frame. For example, if the TF root is camera_link instead of default base_link, then it can be set like this:
<node pkg="rtabmap_ros" type="rgbd_odometry" name="rgbd_odometry" output="screen">


text:
Integrate the function sin^(-1)((2x)/(1+x^2)) | Filo

Class 12

Math

Calculus

Integrals

557
150

Integrate the function

557
150

Connecting you to a tutor in 60 seconds.


text: Mini-Neptunes seem to be common planets. In this work we investigate the possible formation histories and predicted occurrence rates of mini-Neptunes, assuming that the planets form beyond the iceline. We consider pebble and planetesimal accretion accounting for envelope enrichment and two different opacity conditions. We find that the formation of mini-Neptunes is a relatively frequent output when envelope enrichment by volatiles is included, and that there is a “sweet spot” for mini-Neptune formation with a relatively low solid accretion rate of $$\sim 10^{-6}$$ $$M_\oplus$$ $$yr^{-1}$$. This rate is typical for low/intermediate-mass protoplanetary disks and/or disks with low metallicities. With pebble accretion, envelope enrichment and high opacity favor the formation of mini-Neptunes, with more efficient formation at large semimajor axes ($$\gt 30$$ au) and low disk viscosities. For planetesimal accretion, such planets can also form without enrichment, with the opacity being a key aspect in the growth history and favorable formation location. Finally, we show that the formation of Neptune-like planets remains a challenge for planet formation theories.


text:
# Using neural nets to recognize handwritten digits

Neural Networks and Deep Learning

On the exercises and problems

Deep learning

Acknowledgements

Resources

The human visual system is one of the wonders of the world. Consider the following sequence of handwritten digits:

Most people effortlessly recognize those digits as 504192. That ease is deceptive. In each hemisphere of our brain, humans have a primary visual cortex, also known as V1, containing 140 million neurons, with tens of billions of connections between them. And yet human vision involves not just V1, but an entire series of visual cortices - V2, V3, V4, and V5 - doing progressively more complex image processing. We carry in our heads a supercomputer, tuned by evolution over hundreds of millions of years, and superbly adapted to understand the visual world. Recognizing handwritten digits isn't easy. Rather, we humans are stupendously, astoundingly good at making sense of what our eyes show us. But nearly all that work is done unconsciously. And so we don't usually appreciate how tough a problem our visual systems solve.

The difficulty of visual pattern recognition becomes apparent if you attempt to write a computer program to recognize digits like those above. What seems easy when we do it ourselves suddenly becomes extremely difficult. Simple intuitions about how we recognize shapes - "a 9 has a loop at the top, and a vertical stroke in the bottom right" - turn out to be not so simple to express algorithmically. When you try to make such rules precise, you quickly get lost in a morass of exceptions and caveats and special cases. It seems hopeless.

Neural networks approach the problem in a different way. The idea is to take a large number of handwritten digits, known as training examples,

and then develop a system which can learn from those training examples. In other words, the neural network uses the examples to automatically infer rules for recognizing handwritten digits. Furthermore, by increasing the number of training examples, the network can learn more about handwriting, and so improve its accuracy. So while I've shown just 100 training digits above, perhaps we could build a better handwriting recognizer by using thousands or even millions or billions of training examples.

In this chapter we'll write a computer program implementing a neural network that learns to recognize handwritten digits. The program is just 74 lines long, and uses no special neural network libraries. But this short program can recognize digits with an accuracy over 96 percent, without human intervention. Furthermore, in later chapters we'll develop ideas which can improve accuracy to over 99 percent. In fact, the best commercial neural networks are now so good that they are used by banks to process cheques, and by post offices to recognize addresses.

We're focusing on handwriting recognition because it's an excellent prototype problem for learning about neural networks in general. As a prototype it hits a sweet spot: it's challenging - it's no small feat to recognize handwritten digits - but it's not so difficult as to require an extremely complicated solution, or tremendous computational power. Furthermore, it's a great way to develop more advanced techniques, such as deep learning. And so throughout the book we'll return repeatedly to the problem of handwriting recognition. Later in the book, we'll discuss how these ideas may be applied to other problems in computer vision, and also in speech, natural language processing, and other domains.

Of course, if the point of the chapter was only to write a computer program to recognize handwritten digits, then the chapter would be much shorter! But along the way we'll develop many key ideas about neural networks, including two important types of artificial neuron (the perceptron and the sigmoid neuron), and the standard learning algorithm for neural networks, known as stochastic gradient descent. Throughout, I focus on explaining why things are done the way they are, and on building your neural networks intuition. That requires a lengthier discussion than if I just presented the basic mechanics of what's going on, but it's worth it for the deeper understanding you'll attain. Amongst the payoffs, by the end of the chapter we'll be in position to understand what deep learning is, and why it matters.

### Perceptrons

What is a neural network? To get started, I'll explain a type of artificial neuron called a perceptron. Perceptrons were developed in the 1950s and 1960s by the scientist Frank Rosenblatt, inspired by earlier work by Warren McCulloch and Walter Pitts. Today, it's more common to use other models of artificial neurons - in this book, and in much modern work on neural networks, the main neuron model used is one called the sigmoid neuron. We'll get to sigmoid neurons shortly. But to understand why sigmoid neurons are defined the way they are, it's worth taking the time to first understand perceptrons.

So how do perceptrons work? A perceptron takes several binary inputs, x1,x2,, and produces a single binary output:

In the example shown the perceptron has three inputs, x1,x2,x3. In general it could have more or fewer inputs. Rosenblatt proposed a simple rule to compute the output. He introduced weights, w1,w2,, real numbers expressing the importance of the respective inputs to the output. The neuron's output, 0 or 1, is determined by whether the weighted sum jwjxj is less than or greater than some threshold value. Just like the weights, the threshold is a real number which is a parameter of the neuron. To put it in more precise algebraic terms:
output=01if jwjxj thresholdif jwjxj> threshold(1)
That's all there is to how a perceptron works!

That's the basic mathematical model. A way you can think about the perceptron is that it's a device that makes decisions by weighing up evidence. Let me give an example. It's not a very realistic example, but it's easy to understand, and we'll soon get to more realistic examples. Suppose the weekend is coming up, and you've heard that there's going to be a cheese festival in your city. You like cheese, and are trying to decide whether or not to go to the festival. You might make your decision by weighing up three factors:

1. Is the weather good?
2. Does your boyfriend or girlfriend want to accompany you?
3. Is the festival near public transit? (You don't own a car).
We can represent these three factors by corresponding binary variables x1,x2, and x3. For instance, we'd have x1=1 if the weather is good, and x1=0 if the weather is bad. Similarly, x2=1 if your boyfriend or girlfriend wants to go, and x2=0 if not. And similarly again for x3 and public transit.

Now, suppose you absolutely adore cheese, so much so that you're happy to go to the festival even if your boyfriend or girlfriend is uninterested and the festival is hard to get to. But perhaps you really loathe bad weather, and there's no way you'd go to the festival if the weather is bad. You can use perceptrons to model this kind of decision-making. One way to do this is to choose a weight w1=6 for the weather, and w2=2 and w3=2 for the other conditions. The larger value of w1 indicates that the weather matters a lot to you, much more than whether your boyfriend or girlfriend joins you, or the nearness of public transit. Finally, suppose you choose a threshold of 5 for the perceptron. With these choices, the perceptron implements the desired decision-making model, outputting 1 whenever the weather is good, and 0 whenever the weather is bad. It makes no difference to the output whether your boyfriend or girlfriend wants to go, or whether public transit is nearby.

By varying the weights and the threshold, we can get different models of decision-making. For example, suppose we instead chose a threshold of 3. Then the perceptron would decide that you should go to the festival whenever the weather was good or when both the festival was near public transit and your boyfriend or girlfriend was willing to join you. In other words, it'd be a different model of decision-making. Dropping the threshold means you're more willing to go to the festival.

Obviously, the perceptron isn't a complete model of human decision-making! But what the example illustrates is how a perceptron can weigh up different kinds of evidence in order to make decisions. And it should seem plausible that a complex network of perceptrons could make quite subtle decisions:

In this network, the first column of perceptrons - what we'll call the first layer of perceptrons - is making three very simple decisions, by weighing the input evidence. What about the perceptrons in the second layer? Each of those perceptrons is making a decision by weighing up the results from the first layer of decision-making. In this way a perceptron in the second layer can make a decision at a more complex and more abstract level than perceptrons in the first layer. And even more complex decisions can be made by the perceptron in the third layer. In this way, a many-layer network of perceptrons can engage in sophisticated decision making.

Incidentally, when I defined perceptrons I said that a perceptron has just a single output. In the network above the perceptrons look like they have multiple outputs. In fact, they're still single output. The multiple output arrows are merely a useful way of indicating that the output from a perceptron is being used as the input to several other perceptrons. It's less unwieldy than drawing a single output line which then splits.

Let's simplify the way we describe perceptrons. The condition jwjxj>threshold is cumbersome, and we can make two notational changes to simplify it. The first change is to write jwjxj as a dot product, wxjwjxj, where w and x are vectors whose components are the weights and inputs, respectively. The second change is to move the threshold to the other side of the inequality, and to replace it by what's known as the perceptron's bias, bthreshold. Using the bias instead of the threshold, the perceptron rule can be rewritten:

output={01if wx+b0if wx+b>0(2)
You can think of the bias as a measure of how easy it is to get the perceptron to output a 1. Or to put it in more biological terms, the bias is a measure of how easy it is to get the perceptron to fire. For a perceptron with a really big bias, it's extremely easy for the perceptron to output a 1. But if the bias is very negative, then it's difficult for the perceptron to output a 1. Obviously, introducing the bias is only a small change in how we describe perceptrons, but we'll see later that it leads to further notational simplifications. Because of this, in the remainder of the book we won't use the threshold, we'll always use the bias.

I've described perceptrons as a method for weighing evidence to make decisions. Another way perceptrons can be used is to compute the elementary logical functions we usually think of as underlying computation, functions such as AND, OR, and NAND. For example, suppose we have a perceptron with two inputs, each with weight 2, and an overall bias of 3. Here's our perceptron:

Then we see that input 00 produces output 1, since (2)0+(2)0+3=3 is positive. Here, I've introduced the symbol to make the multiplications explicit. Similar calculations show that the inputs 01 and 10 produce output 1. But the input 11 produces output 0, since (2)1+(2)1+3=1 is negative. And so our perceptron implements a NAND gate!

The NAND example shows that we can use perceptrons to compute simple logical functions. In fact, we can use networks of perceptrons to compute any logical function at all. The reason is that the NAND gate is universal for computation, that is, we can build any computation up out of NAND gates. For example, we can use NAND gates to build a circuit which adds two bits, x1 and x2. This requires computing the bitwise sum, x1x2, as well as a carry bit which is set to 1 when both x1 and x2 are 1, i.e., the carry bit is just the bitwise product x1x2:

To get an equivalent network of perceptrons we replace all the NAND gates by perceptrons with two inputs, each with weight 2, and an overall bias of 3. Here's the resulting network. Note that I've moved the perceptron corresponding to the bottom right NAND gate a little, just to make it easier to draw the arrows on the diagram:
One notable aspect of this network of perceptrons is that the output from the leftmost perceptron is used twice as input to the bottommost perceptron. When I defined the perceptron model I didn't say whether this kind of double-output-to-the-same-place was allowed. Actually, it doesn't much matter. If we don't want to allow this kind of thing, then it's possible to simply merge the two lines, into a single connection with a weight of -4 instead of two connections with -2 weights. (If you don't find this obvious, you should stop and prove to yourself that this is equivalent.) With that change, the network looks as follows, with all unmarked weights equal to -2, all biases equal to 3, and a single weight of -4, as marked:
Up to now I've been drawing inputs like x1 and x2 as variables floating to the left of the network of perceptrons. In fact, it's conventional to draw an extra layer of perceptrons - the input layer - to encode the inputs:
This notation for input perceptrons, in which we have an output, but no inputs,
is a shorthand. It doesn't actually mean a perceptron with no inputs. To see this, suppose we did have a perceptron with no inputs. Then the weighted sum jwjxj would always be zero, and so the perceptron would output 1 if b>0, and 0 if b0. That is, the perceptron would simply output a fixed value, not the desired valued (x1, in the example above). It's better to think of the input perceptrons as not really being perceptrons at all, but rather special units which are simply defined to output the desired values, x1,x2,.

The adder example demonstrates how a network of perceptrons can be used to simulate a circuit containing many NAND gates. And because NAND gates are universal for computation, it follows that perceptrons are also universal for computation.

The computational universality of perceptrons is simultaneously reassuring and disappointing. It's reassuring because it tells us that networks of perceptrons can be as powerful as any other computing device. But it's also disappointing, because it makes it seem as though perceptrons are merely a new type of NAND gate. That's hardly big news!

However, the situation is better than this view suggests. It turns out that we can devise learning algorithms which can automatically tune the weights and biases of a network of artificial neurons. This tuning happens in response to external stimuli, without direct intervention by a programmer. These learning algorithms enable us to use artificial neurons in a way which is radically different to conventional logic gates. Instead of explicitly laying out a circuit of NAND and other gates, our neural networks can simply learn to solve problems, sometimes problems where it would be extremely difficult to directly design a conventional circuit.

### Sigmoid neurons

Learning algorithms sound terrific. But how can we devise such algorithms for a neural network? Suppose we have a network of perceptrons that we'd like to use to learn to solve some problem. For example, the inputs to the network might be the raw pixel data from a scanned, handwritten image of a digit. And we'd like the network to learn weights and biases so that the output from the network correctly classifies the digit. To see how learning might work, suppose we make a small change in some weight (or bias) in the network. What we'd like is for this small change in weight to cause only a small corresponding change in the output from the network. As we'll see in a moment, this property will make learning possible. Schematically, here's what we want (obviously this network is too simple to do handwriting recognition!):

If it were true that a small change in a weight (or bias) causes only a small change in output, then we could use this fact to modify the weights and biases to get our network to behave more in the manner we want. For example, suppose the network was mistakenly classifying an image as an "8" when it should be a "9". We could figure out how to make a small change in the weights and biases so the network gets a little closer to classifying the image as a "9". And then we'd repeat this, changing the weights and biases over and over to produce better and better output. The network would be learning.

The problem is that this isn't what happens when our network contains perceptrons. In fact, a small change in the weights or bias of any single perceptron in the network can sometimes cause the output of that perceptron to completely flip, say from 0 to 1. That flip may then cause the behaviour of the rest of the network to completely change in some very complicated way. So while your "9" might now be classified correctly, the behaviour of the network on all the other images is likely to have completely changed in some hard-to-control way. That makes it difficult to see how to gradually modify the weights and biases so that the network gets closer to the desired behaviour. Perhaps there's some clever way of getting around this problem. But it's not immediately obvious how we can get a network of perceptrons to learn.

We can overcome this problem by introducing a new type of artificial neuron called a sigmoid neuron. Sigmoid neurons are similar to perceptrons, but modified so that small changes in their weights and bias cause only a small change in their output. That's the crucial fact which will allow a network of sigmoid neurons to learn.

Okay, let me describe the sigmoid neuron. We'll depict sigmoid neurons in the same way we depicted perceptrons:

Just like a perceptron, the sigmoid neuron has inputs, x1,x2,. But instead of being just 0 or 1, these inputs can also take on any values between 0 and 1. So, for instance, 0.638 is a valid input for a sigmoid neuron. Also just like a perceptron, the sigmoid neuron has weights for each input, w1,w2,, and an overall bias, b. But the output is not 0 or 1. Instead, it's σ(wx+b), where σ is called the sigmoid function* *Incidentally, σ is sometimes called the logistic function, and this new class of neurons called logistic neurons. It's useful to remember this terminology, since these terms are used by many people working with neural nets. However, we'll stick with the sigmoid terminology., and is defined by:
σ(z)11+ez.(3)
To put it all a little more explicitly, the output of a sigmoid neuron with inputs x1,x2,, weights w1,w2,, and bias b is
11+exp(jwjxjb).(4)

At first sight, sigmoid neurons appear very different to perceptrons. The algebraic form of the sigmoid function may seem opaque and forbidding if you're not already familiar with it. In fact, there are many similarities between perceptrons and sigmoid neurons, and the algebraic form of the sigmoid function turns out to be more of a technical detail than a true barrier to understanding.

To understand the similarity to the perceptron model, suppose zwx+b is a large positive number. Then ez0 and so σ(z)1. In other words, when z=wx+b is large and positive, the output from the sigmoid neuron is approximately 1, just as it would have been for a perceptron. Suppose on the other hand that z=wx+b is very negative. Then ez, and σ(z)0. So when z=wx+b is very negative, the behaviour of a sigmoid neuron also closely approximates a perceptron. It's only when wx+b is of modest size that there's much deviation from the perceptron model.

What about the algebraic form of σ? How can we understand that? In fact, the exact form of σ isn't so important - what really matters is the shape of the function when plotted. Here's the shape:

This shape is a smoothed out version of a step function:

If σ had in fact been a step function, then the sigmoid neuron would be a perceptron, since the output would be 1 or 0 depending on whether wx+b was positive or negative* *Actually, when wx+b=0 the perceptron outputs 0, while the step function outputs 1. So, strictly speaking, we'd need to modify the step function at that one point. But you get the idea.. By using the actual σ function we get, as already implied above, a smoothed out perceptron. Indeed, it's the smoothness of the σ function that is the crucial fact, not its detailed form. The smoothness of σ means that small changes Δwj in the weights and Δb in the bias will produce a small change Δoutput in the output from the neuron. In fact, calculus tells us that Δoutput is well approximated by

ΔoutputjoutputwjΔwj+outputbΔb,(5)
where the sum is over all the weights, wj, and output/wj and output/b denote partial derivatives of the output with respect to wj and b, respectively. Don't panic if you're not comfortable with partial derivatives! While the expression above looks complicated, with all the partial derivatives, it's actually saying something very simple (and which is very good news): Δoutput is a linear function of the changes Δwj and Δb in the weights and bias. This linearity makes it easy to choose small changes in the weights and biases to achieve any desired small change in the output. So while sigmoid neurons have much of the same qualitative behaviour as perceptrons, they make it much easier to figure out how changing the weights and biases will change the output.

If it's the shape of σ which really matters, and not its exact form, then why use the particular form used for σ in Equation (3)

σ(z)11+ez
? In fact, later in the book we will occasionally consider neurons where the output is f(wx+b) for some other activation function f(). The main thing that changes when we use a different activation function is that the particular values for the partial derivatives in Equation (5) change. It turns out that when we compute those partial derivatives later, using σ will simplify the algebra, simply because exponentials have lovely properties when differentiated. In any case, σ is commonly-used in work on neural nets, and is the activation function we'll use most often in this book.

How should we interpret the output from a sigmoid neuron? Obviously, one big difference between perceptrons and sigmoid neurons is that sigmoid neurons don't just output 0 or 1. They can have as output any real number between 0 and 1, so values such as 0.173 and 0.689 are legitimate outputs. This can be useful, for example, if we want to use the output value to represent the average intensity of the pixels in an image input to a neural network. But sometimes it can be a nuisance. Suppose we want the output from the network to indicate either "the input image is a 9" or "the input image is not a 9". Obviously, it'd be easiest to do this if the output was a 0 or a 1, as in a perceptron. But in practice we can set up a convention to deal with this, for example, by deciding to interpret any output of at least 0.5 as indicating a "9", and any output less than 0.5 as indicating "not a 9". I'll always explicitly state when we're using such a convention, so it shouldn't cause any confusion.

#### Exercises

• Sigmoid neurons simulating perceptrons, part I
Suppose we take all the weights and biases in a network of perceptrons, and multiply them by a positive constant, c>0. Show that the behaviour of the network doesn't change.

• Sigmoid neurons simulating perceptrons, part II
Suppose we have the same setup as the last problem - a network of perceptrons. Suppose also that the overall input to the network of perceptrons has been chosen. We won't need the actual input value, we just need the input to have been fixed. Suppose the weights and biases are such that wx+b0 for the input x to any particular perceptron in the network. Now replace all the perceptrons in the network by sigmoid neurons, and multiply the weights and biases by a positive constant c>0. Show that in the limit as c the behaviour of this network of sigmoid neurons is exactly the same as the network of perceptrons. How can this fail when wx+b=0 for one of the perceptrons?

### The architecture of neural networks

In the next section I'll introduce a neural network that can do a pretty good job classifying handwritten digits. In preparation for that, it helps to explain some terminology that lets us name different parts of a network. Suppose we have the network:

As mentioned earlier, the leftmost layer in this network is called the input layer, and the neurons within the layer are called input neurons. The rightmost or output layer contains the output neurons, or, as in this case, a single output neuron. The middle layer is called a hidden layer, since the neurons in this layer are neither inputs nor outputs. The term "hidden" perhaps sounds a little mysterious - the first time I heard the term I thought it must have some deep philosophical or mathematical significance - but it really means nothing more than "not an input or an output". The network above has just a single hidden layer, but some networks have multiple hidden layers. For example, the following four-layer network has two hidden layers:
Somewhat confusingly, and for historical reasons, such multiple layer networks are sometimes called multilayer perceptrons or MLPs, despite being made up of sigmoid neurons, not perceptrons. I'm not going to use the MLP terminology in this book, since I think it's confusing, but wanted to warn you of its existence.

The design of the input and output layers in a network is often straightforward. For example, suppose we're trying to determine whether a handwritten image depicts a "9" or not. A natural way to design the network is to encode the intensities of the image pixels into the input neurons. If the image is a 64 by 64 greyscale image, then we'd have 4,096=64×64 input neurons, with the intensities scaled appropriately between 0 and 1. The output layer will contain just a single neuron, with output values of less than 0.5 indicating "input image is not a 9", and values greater than 0.5 indicating "input image is a 9 ".

While the design of the input and output layers of a neural network is often straightforward, there can be quite an art to the design of the hidden layers. In particular, it's not possible to sum up the design process for the hidden layers with a few simple rules of thumb. Instead, neural networks researchers have developed many design heuristics for the hidden layers, which help people get the behaviour they want out of their nets. For example, such heuristics can be used to help determine how to trade off the number of hidden layers against the time required to train the network. We'll meet several such design heuristics later in this book.

Up to now, we've been discussing neural networks where the output from one layer is used as input to the next layer. Such networks are called feedforward neural networks. This means there are no loops in the network - information is always fed forward, never fed back. If we did have loops, we'd end up with situations where the input to the σ function depended on the output. That'd be hard to make sense of, and so we don't allow such loops.

However, there are other models of artificial neural networks in which feedback loops are possible. These models are called recurrent neural networks. The idea in these models is to have neurons which fire for some limited duration of time, before becoming quiescent. That firing can stimulate other neurons, which may fire a little while later, also for a limited duration. That causes still more neurons to fire, and so over time we get a cascade of neurons firing. Loops don't cause problems in such a model, since a neuron's output only affects its input at some later time, not instantaneously.

Recurrent neural nets have been less influential than feedforward networks, in part because the learning algorithms for recurrent nets are (at least to date) less powerful. But recurrent networks are still extremely interesting. They're much closer in spirit to how our brains work than feedforward networks. And it's possible that recurrent networks can solve important problems which can only be solved with great difficulty by feedforward networks. However, to limit our scope, in this book we're going to concentrate on the more widely-used feedforward networks.

### A simple network to classify handwritten digits

Having defined neural networks, let's return to handwriting recognition. We can split the problem of recognizing handwritten digits into two sub-problems. First, we'd like a way of breaking an image containing many digits into a sequence of separate images, each containing a single digit. For example, we'd like to break the image

into six separate images,

We humans solve this segmentation problem with ease, but it's challenging for a computer program to correctly break up the image. Once the image has been segmented, the program then needs to classify each individual digit. So, for instance, we'd like our program to recognize that the first digit above,

is a 5.

We'll focus on writing a program to solve the second problem, that is, classifying individual digits. We do this because it turns out that the segmentation problem is not so difficult to solve, once you have a good way of classifying individual digits. There are many approaches to solving the segmentation problem. One approach is to trial many different ways of segmenting the image, using the individual digit classifier to score each trial segmentation. A trial segmentation gets a high score if the individual digit classifier is confident of its classification in all segments, and a low score if the classifier is having a lot of trouble in one or more segments. The idea is that if the classifier is having trouble somewhere, then it's probably having trouble because the segmentation has been chosen incorrectly. This idea and other variations can be used to solve the segmentation problem quite well. So instead of worrying about segmentation we'll concentrate on developing a neural network which can solve the more interesting and difficult problem, namely, recognizing individual handwritten digits.

To recognize individual digits we will use a three-layer neural network:

The input layer of the network contains neurons encoding the values of the input pixels. As discussed in the next section, our training data for the network will consist of many 28 by 28 pixel images of scanned handwritten digts, and so the input layer contains 784=28×28 neurons. For simplicity I've omitted most of the 784 input neurons in the diagram above. The input pixels are greyscale, with a value of 0.0 representing white, a value of 1.0 representing black, and in between values representing gradually darkening shades of grey.

The second layer of the network is a hidden layer. We denote the number of neurons in this hidden layer by n, and we'll experiment with different values for n. The example shown illustrates a small hidden layer, containing just n=15 neurons.

The output layer of the network contains 10 neurons. If the first neuron fires, i.e., has an output 1, then that will indicate that the network thinks the digit is a 0. If the second neuron fires then that will indicate that the network thinks the digit is a 1. And so on. A little more precisely, we number the output neurons from 0 through 9, and figure out which neuron has the highest activation value. If that neuron is, say, neuron number 6, then our network will guess that the input digit was a 6. And so on for the other output neurons.

You might wonder why we use 10 output neurons. After all, the goal of the network is to tell us which digit (0,1,2,,9) corresponds to the input image. A seemingly natural way of doing that is to use just 4 output neurons, treating each neuron as taking on a binary value, depending on whether the neuron's output is closer to 0 or to 1. Four neurons are enough to encode the answer, since 24=16 is more than the 10 possible values for the input digit. Why should our network use 10 neurons instead? Isn't that inefficient? The ultimate justification is empirical: we can try out both network designs, and it turns out that, for this particular problem, the network with 10 output neurons learns to recognize digits better than the network with 4 output neurons. But that leaves us wondering why using 10 output neurons works better. Is there some heuristic that would tell us in advance that we should use the 10-output encoding instead of the 4-output encoding?

To understand why we do this, it helps to think about what the neural network is doing from first principles. Consider first the case where we use 10 output neurons. Let's concentrate on the first output neuron, the one that's trying to decide whether or not the digit is a 0. It does this by weighing up evidence from the hidden layer of neurons. What are those hidden neurons doing? Well, just suppose for the sake of argument that the first neuron in the hidden layer detects whether or not an image like the following is present:

It can do this by heavily weighting input pixels which overlap with the image, and only lightly weighting the other inputs. In a similar way, let's suppose for the sake of argument that the second, third, and fourth neurons in the hidden layer detect whether or not the following images are present:

As you may have guessed, these four images together make up the 0 image that we saw in the line of digits shown earlier:

So if all four of these hidden neurons are firing then we can conclude that the digit is a 0. Of course, that's not the only sort of evidence we can use to conclude that the image was a 0 - we could legitimately get a


text:
1. Arrange both the polynomial is same order of exponent . It would be good to have terms arrange from highest exponent to lowest exponent i.e
For example
$S(x) = x^2 + 5x^3 + 1 +10x$
$P(x) = x^3 + x^2 + 5 +10x$

Arranging them as per suggestion above

$S(x) = 5x^3 +x^2+ 10x +1$
$P(x) = x^3 + x^2 + 10x +5$

2. Add the like term . By like term ,we mean same exponent terms. There are two method to add like terms. We either can group them horizontally and add it. We can add vertically
Horizontally
$S(x) + P(x)$
$=x^2 + 5x^3 + 1 +10x + (x^3 + x^2 + 5 +10x)$
Arranging both the polynomial in same order
$=5x^3 +x^2+ 10x +1 + (x^3 + x^2 + 10x +5 )$
Opening the parenthesis and grouping the like terms. In addition,no sign are changed when parenthesis are opened
$= 5x^3 + x^3 + x^2+x^2 + 10x +10x + 1+5$
$=6x^3+2x^2+20x+6$

Vertically
$S(x) = x^2 + 5x^3 + 1 +10x$
$P(x) = x^3 + x^2 + 5 +10x$

Arranging both the polynomial in same order

$S(x) = 5x^3 +x^2+ 10x +1$
$P(x) = x^3 + x^2 + 10x +5$

Similarly we can add three or more polynomials

Example-1
$(2x^3 - x+1+x^2) + (x^3 + 6x - 7) + (-3x^2 - 11 + 2x)$
Arranging them in same order
$=(2x^3 + x^2 - x+1) + (x^3 + 6x - 7) + (-3x^2 + 2x - 11)$

Opening the parenthesis and grouping the like terms
$=2x^3+x^3 +x^2-3x^2 -x +6x+2x +1-7-11$
$=3x^3-2x^2+7x-16$

## Subtracting polynomials

We many times need to subtract the two polynomial .It is very similar to adding polynomials only.Subtracting polynomials just means Subtracting the like terms.We need to follow below steps for Subtraction of polynomial
1. Arrange both the polynomial is same order of exponent . It would be good to have terms arrange from highest exponent to lowest exponent i.e
For example
$S(x) = x^2 + 5x^3 + 1 +10x$
$P(x) = x^3 + x^2 + 5 +10x$

Arranging them as per suggestion above

$S(x) = 5x^3 +x^2+ 10x +1$
$P(x) = x^3 + x^2 + 10x +5$

2. Subtract the like term . By like term ,we mean same exponent terms. There are two method to subtract like terms. We either can group them horizontally and subtract it. We can add vertically
Horizontally
$S(x) - P(x)$
$=x^2 + 5x^3 + 1 +10x - (x^3 + x^2 + 5 +10x)$
Arranging both the polynomial in same order
$=5x^3 +x^2+ 10x +1 -(x^3 + x^2 + 10x +5 )$
Opening the parenthesis and grouping the like terms. In Subtract, sign are reversed when parenthesis are opened i.e + becomes - and - becomes +
$= 5x^3 - x^3 + x^2 - x^2 + 10x -10x + 1-5$
Subtract the like term
$=4x^3 - 4$

Vertically
$S(x) = x^2 + 5x^3 + 1 +10x$
$P(x) = x^3 + x^2 + 5 +10x$

Arranging both the polynomial in same order

$S(x) = 5x^3 +x^2+ 10x +1$
$P(x) = x^3 + x^2 + 10x +5$

Similarly we can add three or more polynomials

Example-1
$(2x^3 - x+1+x^2) - (x^3 + 6x - 7) - (-3x^2 - 11 + 2x)$
Arranging them in same order
$=(2x^3 + x^2 - x+1) - (x^3 + 6x - 7) - (-3x^2 + 2x - 11)$

Opening the parenthesis and grouping the like terms
$=2x^3- x^3 +x^2+3x^2 -x -6x-2x +1+7+11$
$=x^3+4x^2-9x+18$


text:
## How should rolling window of parameter estimates look like?

1

I am using Orstein-Uhlenbeck model to model inflation: $$dI_t=\theta(\mu-I_t)dt+\sigma dW_t$$. I have plotted rolling window estimates of all the parameters. However, I do not understand what to conclude from them. I know that this model assumes long-term mean and it would be nice to have convergence is estimates. I would appreciate it if you could help me with conclusions. Red is $$\sigma$$, pink is $$\mu$$ and blue is $$\theta$$.


text:
# Entropy as number of microstates

## Main Question or Discussion Point

As we know S=Q/T. And

Entropy is defined as number of microstates of a system. So does that prove that, the lower the temperature the more the microstates available?

Related Other Physics Topics News on Phys.org
Dale
Mentor
You have to be careful with your equations. $S=k~\ln(\Omega)$ and for a reversible heat transfer $dS=dQ/T$. So at a lower temperature it requires more transfered heat to make the same change in the log of the number of microstates.

Khashishi


text:
# obj-model.h

Go to the documentation of this file.
00001 /*
00002  * obj-model.h
00003  *
00004  * Copyright (C) 2009  Thomas A. Vaughan
00006  *
00007  *
00008  * Redistribution and use in source and binary forms, with or without
00009  * modification, are permitted provided that the following conditions are met:
00010  *     * Redistributions of source code must retain the above copyright
00011  *       notice, this list of conditions and the following disclaimer.
00012  *     * Redistributions in binary form must reproduce the above copyright
00013  *       notice, this list of conditions and the following disclaimer in the
00014  *       documentation and/or other materials provided with the distribution.
00015  *     * Neither the name of the <organization> nor the
00016  *       names of its contributors may be used to endorse or promote products
00017  *       derived from this software without specific prior written permission.
00018  *
00019  * THIS SOFTWARE IS PROVIDED BY THOMAS A. VAUGHAN ''AS IS'' AND ANY
00020  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
00021  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
00022  * DISCLAIMED. IN NO EVENT SHALL THOMAS A. VAUGHAN BE LIABLE FOR ANY
00023  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
00024  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
00025  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
00026  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
00027  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00028  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00029  *
00030  *
00031  * Rendering 3D OBJ models (the Wavefront OBJ format)
00032  */
00033
00034 #ifndef WAVE_GLUT_OBJ_MODEL_H__
00035 #define WAVE_GLUT_OBJ_MODEL_H__
00036
00037 // includes --------------------------------------------------------------------
00038 #include "common/common.h"
00039
00040 #include "glut/renderable.h"
00041 #include "nstream/nstream.h"
00043
00044
00045
00046 namespace obj {
00047
00048
00049 ////////////////////////////////////////////////////////////////////////////////
00050 ///
00051 /// \ingroup models
00052 /// \defgroup obj_model OBJ 3D Models
00053 ///
00054 /// A simple reference implemenation rendering Wavefront OBJ format
00055 /// 3D models.  See http://en.wikipedia.org/wiki/Obj
00056 ///
00057 /// The Open Asset Library implementation also supports OBJ, but there are
00058 /// a few issues there:
00059 ///  - some bugs in how tags are parsed
00060 ///  - only works for native filesystem-based OBJ files
00061 ///
00062 /// OBJ is simple enough so this is a quick reference implementation.
00063 ///
00064 /// A note to users of this model: in the interest of rendering performance,
00065 /// vertex data is copied per vertex.  That is, it is not \em indexed per
00066 /// vertex, but instead the vertex data (including normals etc) is copied.
00067 ///
00068 /// So if you have a large model with a lot of shared vertices, this approach
00069 /// could be too memory intensive.
00070 ///
00071 ////////////////////////////////////////////////////////////////////////////////
00072 /*@{*/
00073
00074
00075
00076 ////////////////////////////////////////////////////////////////////////////////
00077 //
00078 //      public API
00079 //
00080 ////////////////////////////////////////////////////////////////////////////////
00081
00082 class Model : public glut::Renderable {
00083 public:
00084         // virtual destructor --------------------------------------------------
00085         virtual ~Model(void) throw();
00086
00087         // obj::Model class interface methods ----------------------------------
00088
00089         // static factory methods ----------------------------------------------
00090
00091         /// loads an OBJ scene object (aiScene) which can render itself
00092         static smart_ptr<Model> create(IN nstream::Stream * stream,
00093                                 IN float scale);
00094 };
00095
00096
00097
00098 };      // obj namespace
00099
00100
00101 #endif  // WAVE_GLUT_OBJ_MODEL_H__
00102


text:
# some errors in wddm miniport driver

• ### Question

• During my way to build a miniport display driver, I got some errors from the included Dispmprt.h file.

Some of the errors are:

Error C3646 'ACPI_METHOD_ARGUMENT': unknown override specifier MINIPORT C:\Program Files (x86)\Windows Kits\10\Include\10.0.15063.0\shared\acpiioct.h 177
Error    C2065    'ANYSIZE_ARRAY': undeclared identifier  MINIPORT C:\Program Files (x86)\Windows Kits\10\Include\10.0.15063.0\shared\acpiioct.h 169
Error    C3646    'Argument': unknown override specifier MINIPORT C:\Program Files (x86)\Windows Kits\10\Include\10.0.15063.0\shared\acpiioct.h 104

I opened the acpiioct.h and found out that many things like ANYSIZE_ARRAY and _Field_size_bytes_ had not been defined at all.

What can I do to have them solved?

The platform:

- visual studio 2015; windows 7;

I have built other kmdf drivers before with the same configuration and using ntddk.h file and there were no problem.

Monday, December 10, 2018 3:27 PM

### Answers

• I think you are misdiagnosing things here.  ANYSIZE_ARRAY is defined in ntdef.h, which is included by ntddk.h.

However, _Field_size_bytes_ is part of SAL2, the "source annotation language".  These annotations have changed very dramatically from release to release, and it is very easy to get these weird errors if you have a mismatch.  For example, if you are using the latest WDK with a down-level version of the SDK, that will do it.

And, in fact, that's what you have done.  The include file shows you are using the 1703 WDK (build 15063).  That version requires Visual Studio 2017.  You can't use it with Visual Studio 2015.  If you need to use VS2015, then you'll have to fetch the 1607 WDK, build 14393.

Tim Roberts | Driver MVP Emeritus | Providenza & Boekelheide, Inc.

Tuesday, December 11, 2018 1:23 AM