text:
Matching (graph theory)

In the mathematical discipline of graph theory, a matching or independent edge set in an undirected graph is a set of edges without common vertices. Finding a matching in a bipartite graph can be treated as a network flow problem.

Definitions

Given a graph G = (V, E), a matching M in G is a set of pairwise non-adjacent edges, none of which are loops; that is, no two edges share a common vertex.

A vertex is matched (or saturated) if it is an endpoint of one of the edges in the matching. Otherwise the vertex is unmatched.

A maximal matching is a matching M of a graph G that is not a subset of any other matching. A matching M of a graph G is maximal if every edge in G has a non-empty intersection with at least one edge in M. The following figure shows examples of maximal matchings (red) in three graphs.

A maximum matching (also known as maximum-cardinality matching[1]) is a matching that contains the largest possible number of edges. There may be many maximum matchings. The matching number ${\displaystyle \nu (G)}$ of a graph ${\displaystyle G}$ is the size of a maximum matching. Every maximum matching is maximal, but not every maximal matching is a maximum matching. The following figure shows examples of maximum matchings in the same three graphs.

A perfect matching is a matching that matches all vertices of the graph. That is, a matching is perfect if every vertex of the graph is incident to an edge of the matching. Every perfect matching is maximum and hence maximal. In some literature, the term complete matching is used. In the above figure, only part (b) shows a perfect matching. A perfect matching is also a minimum-size edge cover. Thus, the size of a maximum matching is no larger than the size of a minimum edge cover: ν(G) ≤ ρ(G) . A graph can only contain a perfect matching when the graph has an even number of vertices.

A near-perfect matching is one in which exactly one vertex is unmatched. Clearly, a graph can only contain a near-perfect matching when the graph has an odd number of vertices, and near-perfect matchings are maximum matchings. In the above figure, part (c) shows a near-perfect matching. If every vertex is unmatched by some near-perfect matching, then the graph is called factor-critical.

Given a matching M, an alternating path is a path that begins with an unmatched vertex[2] and whose edges belong alternately to the matching and not to the matching. An augmenting path is an alternating path that starts from and ends on free (unmatched) vertices. Berge's lemma states that a matching M is maximum if and only if there is no augmenting path with respect to M.

An induced matching is a matching that is the edge set of an induced subgraph.[3]

Properties

In any graph without isolated vertices, the sum of the matching number and the edge covering number equals the number of vertices.[4] If there is a perfect matching, then both the matching number and the edge cover number are |V | / 2.

If A and B are two maximal matchings, then |A| ≤ 2|B| and |B| ≤ 2|A|. To see this, observe that each edge in B \ A can be adjacent to at most two edges in A \ B because A is a matching; moreover each edge in A \ B is adjacent to an edge in B \ A by maximality of B, hence

${\displaystyle |A\setminus B|\leq 2|B\setminus A|.}$

Further we deduce that

${\displaystyle |A|=|A\cap B|+|A\setminus B|\leq 2|B\cap A|+2|B\setminus A|=2|B|.}$

In particular, this shows that any maximal matching is a 2-approximation of a maximum matching and also a 2-approximation of a minimum maximal matching. This inequality is tight: for example, if G is a path with 3 edges and 4 vertices, the size of a minimum maximal matching is 1 and the size of a maximum matching is 2.

Matching polynomials

A generating function of the number of k-edge matchings in a graph is called a matching polynomial. Let G be a graph and mk be the number of k-edge matchings. One matching polynomial of G is

${\displaystyle \sum _{k\geq 0}m_{k}x^{k}.}$

Another definition gives the matching polynomial as

${\displaystyle \sum _{k\geq 0}(-1)^{k}m_{k}x^{n-2k},}$

where n is the number of vertices in the graph. Each type has its uses; for more information see the article on matching polynomials.

Algorithms and computational complexity

Maximum-cardinality matching

A fundamental problem in combinatorial optimization is finding a maximum matching. This problem has various algorithms for different classes of graphs.

In an unweighted bipartite graph, the optimization problem is to find a maximum cardinality matching. The problem is solved by the Hopcroft-Karp algorithm in time O(VE) time, and there are more efficient randomized algorithms, approximation algorithms, and algorithms for special classes of graphs such as bipartite planar graphs, as described in the main article.

Maximum-weight matching

In a weighted bipartite graph, the optimization problem is to find a maximum-weight matching; a dual problem is to find a minimum-weight matching. This problem is often called maximum weighted bipartite matching, or the assignment problem. The Hungarian algorithm solves the assignment problem and it was one of the beginnings of combinatorial optimization algorithms. It uses a modified shortest path search in the augmenting path algorithm. If the Bellman–Ford algorithm is used for this step, the running time of the Hungarian algorithm becomes ${\displaystyle O(V^{2}E)}$, or the edge cost can be shifted with a potential to achieve ${\displaystyle O(V^{2}\log {V}+VE)}$ running time with the Dijkstra algorithm and Fibonacci heap.[5]

In a non-bipartite weighted graph, the problem of maximum weight matching can be solved in time ${\displaystyle O(V^{2}E)}$ using Edmonds' blossom algorithm.

Maximal matchings

A maximal matching can be found with a simple greedy algorithm. A maximum matching is also a maximal matching, and hence it is possible to find a largest maximal matching in polynomial time. However, no polynomial-time algorithm is known for finding a minimum maximal matching, that is, a maximal matching that contains the smallest possible number of edges.

A maximal matching with k edges is an edge dominating set with k edges. Conversely, if we are given a minimum edge dominating set with k edges, we can construct a maximal matching with k edges in polynomial time. Therefore, the problem of finding a minimum maximal matching is essentially equal to the problem of finding a minimum edge dominating set.[6] Both of these two optimization problems are known to be NP-hard; the decision versions of these problems are classical examples of NP-complete problems.[7] Both problems can be approximated within factor 2 in polynomial time: simply find an arbitrary maximal matching M.[8]

Counting problems

The number of matchings in a graph is known as the Hosoya index of the graph. It is #P-complete to compute this quantity, even for bipartite graphs.[9] It is also #P-complete to count perfect matchings, even in bipartite graphs, because computing the permanent of an arbitrary 0–1 matrix (another #P-complete problem) is the same as computing the number of perfect matchings in the bipartite graph having the given matrix as its biadjacency matrix. However, there exists a fully polynomial time randomized approximation scheme for counting the number of bipartite matchings.[10] A remarkable theorem of Kasteleyn states that the number of perfect matchings in a planar graph can be computed exactly in polynomial time via the FKT algorithm.

The number of perfect matchings in a complete graph Kn (with n even) is given by the double factorial (n − 1)!!.[11] The numbers of matchings in complete graphs, without constraining the matchings to be perfect, are given by the telephone numbers.[12]

Finding all maximally-matchable edges

One of the basic problems in matching theory is to find in a given graph all edges that may be extended to a maximum matching in the graph (such edges are called maximally-matchable edges, or allowed edges). Algorithms for this problem include:

• For general graphs, a deterministic algorithm in time ${\displaystyle O(VE)}$ and a randomized algorithm in time ${\displaystyle {\tilde {O}}(V^{2.376})}$.[13][14]
• For bipartite graphs, if a single maximum matching is found, a deterministic algorithm runs in time ${\displaystyle O(V+E)}$.[15]

Online bipartite matching

The problem of developing an online algorithm for matching was first considered by Richard M. Karp, Umesh Vazirani, and Vijay Vazirani in 1990.[16]

In the online setting, nodes on one side of the bipartite graph arrive one at a time and must either be immediately matched to the other side of the graph or discarded. This is a natural generalization of the secretary problem and has applications to online ad auctions. The best online algorithm, for the unweighted maximization case with a random arrival model, attains a competitive ratio of 0.696.[17]

Characterizations

Kőnig's theorem states that, in bipartite graphs, the maximum matching is equal in size to the minimum vertex cover. Via this result, the minimum vertex cover, maximum independent set, and maximum vertex biclique problems may be solved in polynomial time for bipartite graphs.

Hall's marriage theorem provides a characterization of bipartite graphs which have a perfect matching and the Tutte theorem provides a characterization for arbitrary graphs.

References

1. ^ Alan Gibbons, Algorithmic Graph Theory, Cambridge University Press, 1985, Chapter 5.
2. ^ http://diestel-graph-theory.com/basic.html
3. ^ Cameron, Kathie (1989), "Induced matchings", Special issue for First Montreal Conference on Combinatorics and Computer Science, 1987, Discrete Applied Mathematics, 24 (1–3): 97–102, doi:10.1016/0166-218X(92)90275-F, MR 1011265
4. ^ Gallai, Tibor (1959), "Über extreme Punkt- und Kantenmengen", Ann. Univ. Sci. Budapest. Eötvös Sect. Math., 2: 133–138.
5. ^ Fredman, Michael L.; Tarjan, Robert Endre (1987), "Fibonacci heaps and their uses in improved network optimization algorithms", Journal of the ACM, 34 (3): 596–615, doi:10.1145/28869.28874
6. ^ Yannakakis, Mihalis; Gavril, Fanica (1980), "Edge dominating sets in graphs" (PDF), SIAM Journal on Applied Mathematics, 38 (3): 364–372, doi:10.1137/0138030.
7. ^ Garey, Michael R.; Johnson, David S. (1979), Computers and Intractability: A Guide to the Theory of NP-Completeness, W.H. Freeman, ISBN 0-7167-1045-5. Edge dominating set (decision version) is discussed under the dominating set problem, which is the problem GT2 in Appendix A1.1. Minimum maximal matching (decision version) is the problem GT10 in Appendix A1.1.
8. ^ Ausiello, Giorgio; Crescenzi, Pierluigi; Gambosi, Giorgio; Kann, Viggo; Marchetti-Spaccamela, Alberto; Protasi, Marco (2003), Complexity and Approximation: Combinatorial Optimization Problems and Their Approximability Properties, Springer. Minimum edge dominating set (optimisation version) is the problem GT3 in Appendix B (page 370). Minimum maximal matching (optimisation version) is the problem GT10 in Appendix B (page 374). See also Minimum Edge Dominating Set and Minimum Maximal Matching in the web compendium.
9. ^ Leslie Valiant, The Complexity of Enumeration and Reliability Problems, SIAM J. Comput., 8(3), 410–421
10. ^ Bezáková, Ivona; Štefankovič, Daniel; Vazirani, Vijay V.; Vigoda, Eric (2008). "Accelerating Simulated Annealing for the Permanent and Combinatorial Counting Problems". SIAM Journal on Computing. 37 (5): 1429–1454. CiteSeerX 10.1.1.80.687. doi:10.1137/050644033.
11. ^ Callan, David (2009), A combinatorial survey of identities for the double factorial, arXiv:0906.1317, Bibcode:2009arXiv0906.1317C.
12. ^ Tichy, Robert F.; Wagner, Stephan (2005), "Extremal problems for topological indices in combinatorial chemistry" (PDF), Journal of Computational Biology, 12 (7): 1004–1013, doi:10.1089/cmb.2005.12.1004, PMID 16201918.
13. ^ Rabin, Michael O.; Vazirani, Vijay V. (1989), "Maximum matchings in general graphs through randomization", Journal of Algorithms, 10 (4): 557–567, doi:10.1016/0196-6774(89)90005-9
14. ^ Cheriyan, Joseph (1997), "Randomized ${\displaystyle {\widetilde {O}}(M(|V|))}$ algorithms for problems in matching theory", SIAM Journal on Computing, 26 (6): 1635–1655, doi:10.1137/S0097539793256223
15. ^ Tassa, Tamir (2012), "Finding all maximally-matchable edges in a bipartite graph", Theoretical Computer Science, 423: 50–58, doi:10.1016/j.tcs.2011.12.071
16. ^ Karp, Richard M.; Vazirani, Umesh V.; Vazirani, Vijay V. (1990). "An optimal algorithm for on-line bipartite matching" (PDF). Proceedings of the 22nd Annual ACM Symposium on Theory of Computing (STOC 1990). pp. 352–358. doi:10.1145/100216.100262.
17. ^ Mahdian, Mohammad; Yan, Qiqi (2011). "Online bipartite matching with random arrivals: an approach based on strongly factor-revealing LPs". Proceedings of the Forty-Third Annual ACM Symposium on Theory of Computing. pp. 597–606. doi:10.1145/1993636.1993716.
18. ^ See, e.g., Trinajstić, Nenad; Klein, Douglas J.; Randić, Milan (1986), "On some solved and unsolved problems of chemical graph theory", International Journal of Quantum Chemistry, 30 (S20): 699–742, doi:10.1002/qua.560300762.


text:
# intersection_plane_plane

compas.geometry.intersection_plane_plane(plane1, plane2, tol=1e-06)[source]

Computes the intersection of two planes

Parameters
• plane1 (tuple) – The base point and normal (normalized) defining the 1st plane.

• plane2 (tuple) – The base point and normal (normalized) defining the 2nd plane.

• tol (float, optional) – A tolerance for membership verification. Default is 1e-6.

Returns

line (tuple) – Two points defining the intersection line. None if planes are parallel.

format: plain


text:
US Medicare Payment Data - comparison of payments among states¶

In this script we compare Medicare payments in 2012 between several states. There are many ways to do this. For example, we could look at the total or average of all payments to all providers in a state. However the total payment at the state level would scale with the size of the state, and the average payment may be affected by the composition of patients and providers within the state. To reduce these effects, we can look at the average payment made within a state to all providers of a specific type, or for all services of a particular type. These averages can then be compared between states.

The CMS payment data can be obtained here.

Here is some documentation about the data set.

First we import the libraries that we will use in the script.

In [1]:
%matplotlib inline

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

Next we read in the data for each state that we intend to analyze. We use a dictionary comprehension to create a dictionary that maps the two-letter state name string to the Pandas DataFrame object for the state.

In [2]:
state_names = ["FL", "MI"]
fnames = [x + "-subset.csv.gz" for x in state_names]
state_data = {k: pd.read_csv(f, compression="gzip") for (k, f) in zip(state_names, fnames)}
print([type(state_data[k]) for k in state_names])
print([state_data[k].shape for k in state_names])

[<class 'pandas.core.frame.DataFrame'>, <class 'pandas.core.frame.DataFrame'>]
[(667995, 27), (339456, 27)]

We plan to focus on individual providers here, so we exclude the provider organizations.

In [3]:
for k in state_names:
ix = (state_data[k].loc[:, "nppes_entity_code"] == "I")
state_data[k] = state_data[k].loc[ix, :]
print(state_data[k].shape)

(636249, 27)
(329457, 27)

Comparing mean payments by provider type¶

We start by calculating the total medicare payment for each type of service within each provider, within each state.

In [4]:
for k in state_names:
state_data[k]["total_payment"] = state_data[k]["line_srvc_cnt"] * state_data[k]["average_Medicare_payment_amt"]
print(state_data[k].shape)

(636249, 28)
(329457, 28)

The simplest comparison is to look at the total of all payments by state. We calculate this quantity in the next cell (in billions of US dollars). If we are comparing Florida to Michigan, Florida's total population is 19.3 million and Michigan's is 9.8 million. Florida's population of people over 65 years old is 4.1 million, and Michigan's is 1.5 million.

In [5]:
tot = {k : state_data[k].loc[:, "total_payment"].sum() / 1e9 for k in state_names}
print(tot)
print("FL/MI payent ratio: %.3f" % (tot["FL"] / tot["MI"]))
print("FL/MI population ratio: %.3f" % (19.3 / 9.8))
print("Fl/MI 65+ population ratio: %.3f" % (4.1 / 1.5))

{'FL': 6.5373885239300451, 'MI': 2.433227685210483}
FL/MI payent ratio: 2.687
FL/MI population ratio: 1.969
Fl/MI 65+ population ratio: 2.733

Next we aggregate the data by provider, taking the sum of all payments to each provider within each state. Again we use dictionary comprehensions, so the result is a dictionary mapping state abbreviations to a DataFrames containing the aggregated values. We aggregate the payments by summing, and we aggregate the provider type by taking the first value within the group (since it is constant within each group). The first function extracts the first element from a Pandas Series object.

In [6]:
first = lambda x: x.iloc[0]

by_provider = {k: state_data[k].groupby("npi").agg({"total_payment": np.sum, "provider_type": first}) for k in state_names}
print({k: type(by_provider) for k in state_names})
print({k: by_provider[k].shape for k in state_names})

{'FL': <class 'dict'>, 'MI': <class 'dict'>}
{'FL': (49076, 2), 'MI': (31138, 2)}

This is what we have at this point:

In [7]:
print(state_names[0])

FL
total_payment       provider_type
npi
1003000381       16411.12  Physical Therapist
1003000522       97836.25     Family Practice
1003000720        7019.95  Nurse Practitioner
1003002817      129028.63     Family Practice

Next we aggregate again, this time aggregating over all providers of a given type.

In [8]:
by_provider_type = {k: by_provider[k].groupby("provider_type").agg(np.mean) for k in state_names}
print({k: type(by_provider_type[k]) for k in state_names})
print({k: by_provider_type[k].shape for k in state_names})

{'FL': <class 'pandas.core.frame.DataFrame'>, 'MI': <class 'pandas.core.frame.DataFrame'>}
{'FL': (75, 1), 'MI': (73, 1)}

A few provider types are used in one state but not the other:

In [9]:
ptypes = {k: by_provider_type[k].index for k in state_names}
print(set(ptypes["MI"]) - set(ptypes["FL"]))
print(set(ptypes["FL"]) - set(ptypes["MI"]))

{'Unknown Supplier/Provider'}
{'Unknown Physician Specialty Code', 'Independent Diagnostic Testing Facility', 'Peripheral Vascular Disease'}

This is what these data look like for one state, the numbers are the average total payment to a provider of a given type within the state:

In [10]:
by_provider_type["MI"].head()

Out[10]:
total_payment
provider_type
Allergy/Immunology 53693.749917
Anesthesiologist Assistants 3732.635000
Anesthesiology 36701.068349
Audiologist (billing independently) 5019.827558

Up to this point, we had a separate data set for the each state. Now that we have reduced the data to the main quantity of interest, we can merge the state-level data sets into a single data set, to make it easier to compare the states.

In [11]:
for k in state_names:
by_provider_type[k].columns = [k]
by_provider_type = pd.concat(list(by_provider_type.values()), axis=1)

Here is a peek at the merged data set:

In [12]:
print(by_provider_type.head())

                                                FL             MI
Allergy/Immunology                   102898.723737   53693.749917
Anesthesiologist Assistants           11456.243731    3732.635000
Anesthesiology                        41096.014748   36701.068349
Audiologist (billing independently)   12976.023042    5019.827558

Since some provider types are not present in a given state, after merging their payment value will be nan. We can confirm that here:

In [13]:
print(pd.isnull(by_provider_type).sum())

FL    1
MI    3
dtype: int64

In this setting, it might be reasonable to replace the missing values with zero:

In [14]:
by_provider_type = by_provider_type.replace(np.nan, 0)
print(pd.isnull(by_provider_type).sum())

FL    0
MI    0
dtype: int64

Next we make a scatterplot showing the average payments by provider type in one state versus the other.

In [15]:
plt.plot(by_provider_type.loc[:, state_names[0]], by_provider_type.loc[:, state_names[1]], 'o', alpha=0.6)
plt.plot([0, 1e6], [0, 1e6], '-', color='grey')
plt.xlabel(state_names[0] + " provider types", size=15)
plt.ylabel(state_names[1] + " provider types", size=15)

Out[15]:
<matplotlib.text.Text at 0x7fc871ba74a8>

Since the payment levels are quite skewed, we will plot them on a log scale. Note that the zeros that were introduced into the dataset will become nan after taking the log, and hence will be dropped from the plot. This is OK because those probably reflect administrative differences in how the coding is handled, rather than being true zeros.

In [16]:
x = np.log(by_provider_type.loc[:, state_names[0]]) / np.log(2)
y = np.log(by_provider_type.loc[:, state_names[1]]) / np.log(2)
plt.plot(x, y, 'o', alpha=0.6)
plt.plot([10, 22], [10, 22], '-', color='grey')
plt.xlabel(state_names[0] + " provider types (log2)", size=15)
plt.ylabel(state_names[1] + " provider types (log2)", size=15)

Out[16]:
<matplotlib.text.Text at 0x7fc86ffa4780>

We can also drill down and examine some of the provider types with the greatest difference in payment levels between the two states.

In [17]:
ii = (y - x > 2)
print(by_provider_type.loc[ii, :])

                                     FL           MI
Unknown Supplier/Provider      0.000000  115743.6225

In [18]:
ii = (x - y > 2)
print(by_provider_type.loc[ii, :])

                                                    FL         MI
Independent Diagnostic Testing Facility   15037.770000      0.000
Nuclear Medicine                         281914.943600  70312.875
Peripheral Vascular Disease              126513.343333      0.000
Unknown Physician Specialty Code          90887.808889      0.000

Comparison by service type¶

In [19]:
by_service_type = {k: state_data[k].groupby("hcpcs_code").agg({"average_Medicare_payment_amt": np.mean, "hcpcs_description": first}) for k in state_names}

            average_Medicare_payment_amt         hcpcs_description
hcpcs_code
00103                         136.483071     Anesth blepharoplasty
00104                          64.475266       Anesth electroshock
00120                         142.361123        Anesth ear surgery
00126                          60.844473        Anesth tympanotomy
00140                         128.521658  Anesth procedures on eye

Next we concatenate the state-level results into a single data frame. First we rename the variable of interest to the state's name.

In [20]:
for k in state_names:
by_service_type[k].rename(columns = {"average_Medicare_payment_amt": k}, inplace=True)
by_service_type = pd.concat(list(by_service_type.values()), axis=1)
print(by_service_type.shape)

(3903, 4)
FL         hcpcs_description          MI  \
00103  136.483071     Anesth blepharoplasty   93.955017
00104   64.475266       Anesth electroshock   46.538973
00120  142.361123        Anesth ear surgery  150.657667
00126   60.844473        Anesth tympanotomy   44.884947
00140  128.521658  Anesth procedures on eye  100.105619

hcpcs_description
00103     Anesth blepharoplasty
00104       Anesth electroshock
00120        Anesth ear surgery
00126        Anesth tympanotomy
00140  Anesth procedures on eye

After concatenating, the hcpcs_description column will be duplicated, since this column was present in the data file for each state. We don't really need to remove the extra copy, but we can.

In [21]:
ii = np.flatnonzero(by_service_type.columns == "hcpcs_description")
jj = [0, 1, 2, 3]
jj.remove(ii[1])
by_service_type = by_service_type.iloc[:, jj]

               FL         hcpcs_description          MI
00103  136.483071     Anesth blepharoplasty   93.955017
00104   64.475266       Anesth electroshock   46.538973
00120  142.361123        Anesth ear surgery  150.657667
00126   60.844473        Anesth tympanotomy   44.884947
00140  128.521658  Anesth procedures on eye  100.105619

In [22]:
plt.plot(by_service_type.loc[:, state_names[0]], by_service_type.loc[:, state_names[1]], 'o', alpha=0.6)
plt.plot([0, 1e4], [0, 1e4], '-', color='grey', alpha=0.5)
plt.xlabel(state_names[0] + " service types", size=15)
_ = plt.ylabel(state_names[1] + " service types", size=15)

In [23]:
x = np.log(by_service_type.loc[:, state_names[0]]) / np.log(2)
y = np.log(by_service_type.loc[:, state_names[1]]) / np.log(2)
plt.plot(x, y, 'o', alpha=0.6)
plt.plot([0, 14], [0, 14], '-', color='grey', alpha=0.5)
plt.xlim(0, 14)
plt.ylim(0, 14)
plt.xlabel(state_names[0] + " service types", size=15)
_ = plt.ylabel(state_names[1] + " service types", size=15)

Another way to compare the two states is to examine the distribution of provider total payments. A quantile-quantile (QQ) plot shows in detail how two distributions differ.

In [24]:
provider_totals = {k: state_data[k].groupby("npi").agg({"total_payment": np.sum}) for k in state_names}

            total_payment
npi
1003000605       19619.11
1003001884       28618.18
1003003021      107302.65
1003005315       67463.36
1003006503       35688.16

We will need the a set of percentiles from the two data sets.

In [25]:
log_provider_totals = {k : np.log10(provider_totals[k]) for k in state_names}
q = np.linspace(1, 99, 50).tolist()
pctls = {k: np.percentile(log_provider_totals[k], q) for k in state_names}
q = np.linspace(10, 91, 9).tolist()
dcls = {k: np.percentile(log_provider_totals[k], q) for k in state_names}

Here is the QQ plot. It shows that Michigan providers receive lower payments than Florida providers at every point in the distribution.

In [26]:
plt.plot(pctls[state_names[0]], pctls[state_names[1]], '-', lw=4)
plt.plot(dcls[state_names[0]], dcls[state_names[1]], 'o', ms=10, mec='blue', mfc='white')
plt.plot([2, 7], [2, 7], '-', color='grey')
plt.grid(True)
plt.xlabel("log10 " + state_names[0] + " payment quantiles", size=15)
_ = plt.ylabel("log10 " + state_names[1] + " payment quantiles", size=15)

Exercises:

• Compare the states' total Medicare payments in terms of the number of unique beneficiaries, and the number of unique beneficiary days (see the data documentation for the variable names). Interpret the results relative to the states' population sizes and Medicare eligible (65 years of age or older) population sizes.

• Calculate the proportion of providers in each state who belong to each provider type. Compare these proportions between the two states using some of the techniques demonstrated above.

• Repeat some of the analyses above using alternative measures to count the services (number of unique beneficiaries and number of unique beneficiary days).


text:
# Integral for Bhattacharyya distance between two Cauchy distributions

I need to perform the following integral to calculate the Bhattacharyya distance between two Cauchy distributions:

$$I = \frac{\sqrt{b_+ b_-}}{\pi}\int_{-\infty}^{\infty}dx\,\frac{1}{\sqrt{\left[(x-1)^2+b_+^2\right]\left[(x+1)^2+b_-^2\right]}}$$

Here, $$b_+ >0$$, $$b_->0$$, and $$b_+ \neq b_-$$. The case $$b_+ = b_-$$ is solved here.

I saw on another forum that Maple can apparently express this integral in terms for Elliptic functions. But the user displayed their answer in a form that is completely unreadable. I don't have access to Maple, so I tried reproducing the result with Mathematica 12.1 but it seems incapable to find an answer.

Anyone know how I can help Mathematica do this? The code that is not sufficient is:

int = (Sqrt[bp bm]/Pi) Integrate[1/Sqrt[((x - 1)^2 + bp^2) ((x + 1)^2 + bm^2)], {x, -Infinity, Infinity}, Assumptions -> {bp > 0, bm > 0}]

• With Maple I have: (2 Sqrt[2] Sqrt[bm bp] EllipticK[(2 (bm^4 - 2 bm^2 (-4 + bp^2) + (4 + bp^2)^2))/( bm^4 - 2 bm^2 (-4 + bp^2) + (4 + bp^2)^2 + (4 + bm^2 + bp^2) Sqrt[ bm^4 - 2 bm^2 (-4 + bp^2) + (4 + bp^2)^2])])/(Sqrt[ 4 + bm^2 + bp^2 + Sqrt[ bm^4 - 2 bm^2 (-4 + bp^2) + (4 + bp^2)^2]] \[Pi]) – Mariusz Iwaniuk Dec 17 '20 at 14:24
• I confirm that it works numerically. Thanks a lot. – Ben Dec 17 '20 at 14:34

$$I = \frac{4}{\pi}\sqrt{\frac{b_+ b_-}{A_+}} K\left(\frac{A_-}{A_+}\right).$$

Here, I have introduced:

$$A_\pm = 4 + (b_+ \pm b_-)^2,$$

and $$K$$ is the complete elliptic integral of the first kind as defined by Mathematica's EllipticK[].

Edit: I have incorporated Andrea's further simplifications.

• you could complete simplification with ((4*Sqrt[bm*bp])/(Sqrt[ap]*Pi))*EllipticK[am/ap] – Andreas Dec 17 '20 at 18:15
• I checked it numerically and I don't think it gives the right result? There might be other simplifications though. – Ben Dec 17 '20 at 19:32
• {((4*Sqrt[bp*bm])/(Pi*(Sqrt[ap] - Sqrt[am])))* EllipticK[4/(2 - Sqrt[ap/am] - Sqrt[am/ap])] /. {ap -> 4 + (bp + bm)^2, am -> 4 + (bp - bm)^2}, ((4*Sqrt[bm*bp])/(Sqrt[ap]*Pi))* EllipticK[am/ap] /. {ap -> 4 + (bp + bm)^2, am -> 4 + (bp - bm)^2}} /. {bp -> 0.7, bm -> 0.3} – Andreas Dec 17 '20 at 19:41
• You're correct, it works. I had made a mistake. I will change it. – Ben Dec 17 '20 at 20:11


text:
## Saturday, 20 April 2013

### Check out my quora blog, too

To see newer posts, check out inference.quora.com

## Friday, 5 October 2012

### The big data business model

This is of course in reference to this clip from South Park:

## Monday, 18 June 2012

### How to think outside the box - literally

This blog post is about innovation and creativity. More precisely it's about a stupid little technique, that once helped me be creative and achieve something unexpected: come up with a quantum physics technique, and publish a paper, with barely any training in physics.

Most big ideas and breakthroughs are unexpected. They are not part of your plan for next week, in case of a company, probably not on your roadmap for next year. It is important to have a roadmap, but it is just as crucial to realise that your potentially biggest achievements are not yet on it. Keeping that in mind, what can you do right now to trigger that next big idea coming? This may be your moment: follow these instructions. It may take some time, but bear with me. I mean, really, try this, if you can't commit yourself to try something creative right now, when you have nothing better to do than read this crappy blog post, then when will you?

 My thinking box
Step 1: To be creative, you want to go beyond your confort zone: you want to think outside the box. What's the first thing you need for that? A box. That's right. If you want to think outside the box, you have to have the box first. So go and get one now. A cardboard box, large or small, will easily do. I used the one pictured on the right, which I found in the recycling bin in the office. Resume reading when you found your box.

Step 2: Now that you have the box, what's inside? Inside the box are things that you are confortable with thinking about. If you have a to do list, a project plan, a roadmap for next year(s), those all go inside the box. Everything that you anticipate to happen to you and your work, everything predictable, those all go inside the box.

Step 3+: When you're ready and have the box, you can start thinking outside of it. You will likely need a pencil, to take notes, I simply sketch stuff on the wall of the box. You may prefer a quiet place to think: this time, try a busy place: remember, you're outside your comfort zone now. The rest is up to you. Let me know if you how you saved the world using this technique.

## Wednesday, 6 June 2012

### My big fat data-driven wedding: how to decide when to get married

There's all this hype about the importance of using data to support critical business decisions. But few things can be more important in someone's life than choosing the partner they plan to spend the rest of their lives with. So look no further for a data-driven approach to support this critical decision: when should you stop looking and finally get married?The formal setup of the problem is as follows: there are $$N$$ potential marriage partners out there, and I want to choose the one and only partner of my life. My goal of course is to settle with no-one but the very best partner to marry.I can 'date' each candidate sequentially in a random order. When I date someone, I can evaluate how she ranks in comparison to all previous partners I dated. Based on this, I can either go on and marry her in which case the process stops, or reject her and carry on dating more partners. Once a partner was rejected she cannot be called back again. The problem is hard, because it's assumed that I know nothing about a partner before I actually dated her. Some might immediately recognise the typical problem of exploration versus expoitation hidded in here. So when should I stop dating and marry someone?Let's say I want to maximise the probability of ending up with the best partner as wife. Firstly, there's no point in selecting someone who's not the best so far, because that also implies she cannot be the best overall. Also, the earlier you marry, the higher the chance of missing out on someone without even dating her. It turns out the optimal strategy is as follows: date the first $$k$$ candidates and reject them irrespective of their 'performance'. Then, keep on dating partners, and choose the first one, who ranks above everyone I dated before. The threshold $$k$$ depends on $$N$$, assymptotically it's about a third of all candidates ( $$k \approx \frac{1}{e}N$$ see slides for details).I gave a talk about this theory in our lab at Cambridge not long after getting married, and did a quick back-of-the-envelope calculation to see how close my marriage decision was to the optimal. The theory predicts based on my criteria I should've dated roughly 196,132 partners before getting married. Well, I leave it to you to guess whether I did that or not :) But at the end of the day I ended up making the best decision, so it doesn't matter.I got recently reminded of this problem, as I was thinking about making instant decisions in the PeerPerks platform, which as it turns out  can be modelled as a more complicated version of the optimal marriage problem, called probabilistic submodular multiple-choice secretary problem. In PeerPerks, we allow brands to offer freebies to influential users of social media. Our job is to make sure the promotion generates maximal success based on various indicators. In PeerPerks, every once in a while a user turns up on the website and authenticates using their twitter or facebook. At this point we have access to their social media profile, so it's like we've dated that person: we have a little insight into which topics she's interested in, what kind of people she's connected to, etc. Based on this instant evaluation our algorithm decides whether or not the person qualifies to receive a perk or not.

## Tuesday, 6 March 2012

### Sock gnomes versus laws of large numbers

Whilst finishing one of my previous posts on causality and influence, I was pairing up socks that just came out of the laundry. I was outraged to observe that, yet again, I was hardly able to pair up over half of the socks. 'This cannot happen by chance, we must have sock gnomes, who deliberately mess up your socks' - I was thinking to myself. Clearly, I'm not alone with this problem:

But then I realised I may be wrong, maybe it CAN happen by chance. Maybe it's in fact EXPECTED to happen, and if I could perfectly pair up all socks, I ought to be more surprised. So I started formulating the problem:
We have $$n$$ perfectly distinguishable pairs of socks in our laundry bag. We decide to do our laundry, but because our washing machine has finite capacity, we can only wash $$k \leq 2n$$ pieces of socks. We select the particular socks to be washed by drawing from the laundry bag uniformly at random, without replacement. After this procedure in our sample of $$k$$ socks there will be a number $$Y_{n,k}$$ of them that won't have their pairs in the sample. What is the expected value of  $$Y_{n,k}$$ as a function of $$n$$ and $$k$$?
First I wanted to give an analytic answer, but I had waaay to many other things to use my brainpower for, so sorry, this blog post is based on Monte Carlo estimates, that I could quickly hack together. (For those of you who doesn't know, Monte Carlo is a fancy name for simulating a bunch of random outcomes and then taking averages)

I came up with a visualisation to explain what's going on. Here's a simple (and ugly) version of plot if we have only $$n=2$$ pairs of socks, with explanation below.

2 pairs of socks means a total of 4 pieces in our laundry bag that we may select to wash. The x axis shows the number of socks washed (k), which therefore ranges between 0 and 4. The y axis shows the number of socks left unpaired. The gray rectangles show the distribution of unpaired socks with white meaning probability 0, black meaning probability 1 (see colorbar on the right).
If we wash 0 socks, with probability 1 we have 0 unpaired socks, hence the black rectangle at (0,0), meaning if we wash 0 socks, with probability 1 we have 0 unpaired socks.
If we wash a single sock, we obviously can't pair it up with anything, thus with probability one, we have 1 unpaired sock, hence the black rectangle at (1,1).
If we wash 2 pieces, than with probability ~0.65 those two won't be pairs, so we're left with two unpaired socks, hence the dark gray rectangle at (2,2); with probability ~0.35 we can actually pair them up and we're left with 0 unpaired socks, hence the lighter gray rectangle at (2,0). It's impossible to have a single unpaired sock, hence (2,1) is white.
And I hope you get it by now and could finish the reasoning all the way to the end, and figure out why the graph is always symmetric.

The blue dotted line shows the maximum of the number of unpaired socks, it's easy to show you can never have more than $$min(k,2n-k)$$ unpaired socks. So everything above the blue dotted line should be white. The red line shows the expected number of unpaired socks as a function of $$k$$, the number of pieces we have washed (remember, $$n$$ is now fixed)

Below is the same kind of graph for the case we have $$n=10$$ pairs of socks. Notice the checkerboard-like pattern which due to the fact that the parity of the number of unpaired socks always corresponds to the parity of the number of socks washed. Pretty pleasing to the eye, isn't it?
The expected number of socks starts to be a smoother function of the number of socks washed, and it becomes even further as we further increase this number to 50, as the following figure shows:
The shape of the expectation looks the same as a function of $$\frac{k}{2n}$$, the variance around this mean however has decreased substantially.
If we further increase the number of socks to 200 pairs, the variance practically shrinks to 0:
And for a little animation putting all this together for growing values of $$n$$ click here.

So what's the answer? Well, based on Monte Carlo simulations I conjecture that the following law of large numbers holds for sock pairing:

As $$n \rightarrow \infty$$ and $$k \rightarrow \infty$$ but such that $$\lambda = \frac{k}{n2}$$ is fixed, the fraction of unpaired socks converges to a deterministic value with probability one:

$\frac{Y_{n,k}}{2n} \stackrel{1}{\rightarrow} \lambda\left(1-\lambda\right)$

So, the expected fraction of socks that you can't pair up can be as high as half of the total number of socks you have washed, so my analysis allowed me to rule out the alternative hypothesis of sock gnomes messing up my socks.

Take home message #1: Before you start believing in sock gnomes, satan, god, and other supernatural powers, check if events that seem unlikely really cannot be explained by more parsimonious models of probability.

Take home message #2: If you want to be able to pair up your socks, buy loads of the same colour and you won't have such problems.

## Sunday, 29 January 2012

### Flattr: should we use it in academia?

I have just found the service that I have been looking for for a while: Flattr. It's a social micropayment platform, an evolution of "Buy me a coffee" buttons that started appearing on blogs a few years ago.
It allows you to donate and distribute a small amount of money each month to authors of content that you "liked" that month.

Here's how it works: you sign up on Flattr, top up your account with a small amount, let's say $5 a month, then you come back to my blog and press the Flattr button on the right. At the end of the month, your$5 will be distributed equally to authors of things you Flattr'ed, so I get my share of your \$5. This way, the small contributions add up, and the community can reward authors of content they frequently read and who inspire them. I personally think, details and implementation apart, this general idea is brilliant. I'm definitely going to top up my account and suggest people, whose stuff I read to add a flattr button.

I actually first thought of such a social redistribution system in terms of academic funding, and almost launched a buy me a coffee button on my publications page to make a point. Now, having found Flattr, I started thinking again about how could such a system be used for a more equal redistribution  of funding in academia. Let's review some current inefficiencies in academic research funding:

Inflexible fellowship values: Money in academia is often awarded in discrete chunks, especially for junior positions. You apply for various grants, fellowships, scholarships, and you either get the whole thing for x years or you get nothing. If you for example got a Gates Scholarship to study for a PhD in Cambridge, you got it for three years, the amount you get is practically fixed for 3 years (other than following inflation), you can't cease to be a Gates scholar if you do poorly, you can't become one if you perform better. It's pretty much fixed and independent of your performance for 3-4 years. This lack of adaptation may generate random inequalities: there are some fellowships that are only awarded every second year for example; it's close to impossible to get a decent academic salary in certain countries, independently of the quality of research being done.

Poor metrics: In my opinion academia struggles with a big unsolved problem: How to attribute authority, how to measure and reward someone's contribution, inspiration, impact? The most widespread metrics are number of citations, h-index, impact factors,
eigenfactorsetc. These are crude, heuristic measures of influence and authority. They certainly make some intuitive sense and vaguely correlate with what we really wanted to measure, but these are not a solution, just complicated features. If we want a fairer system we have to establish new standards for measuring influence/importance/etc. (Well, first of all we have to agree what we want to measure)

Discrete random variables: When measuring output, there are way too many discrete variables involved: Is your paper accepted? Yes or No? Does this paper cite relevant papers of yours appropriately? Yes or No? Which journal did your stuff appear in? These binary decisions rectify the opinion formed about your paper: if your paper is accepted, it means some undisclosed set of reviewers who looked at it didn't find a mistake and liked it. But how much did they like it? Were they absolutely amazed by it, or did they just not care enough to reject it? Were the reviewers respected academics, or were they incompetent? What does the rest of the community think? If I want to state that in my opinion a particular piece of research represents huge step forward in machine learning, pretty much my only option to measurably reward the authors is to cite their paper, but maybe my work is on something completely different, so it's inappropriate.

All in all, the measurement of scientific output is inherently inefficient because of these binary/discrete decisions involved. Everybody with engineering background knows that discretisation means loss of information: measurement of academic output is discretised rather arbitrarily, with little consideration given to effective information transfer. And this loss of bandwith constrains the system of academic rewards to be either slow and accurate (you have to grow old to get recognition) or fast but inaccurate (you will be judged on the basis of a few noisy random variables). Often I think it's both slow and inaccurate :)

Does Flattr suggest a solution? Can something like Flattr help improve academic funding and assignment of credits? Imagine you get a fellowship, and you have to re-allocate 5% of your earnings each month to fellow scientists according to how important their contributions are to your field? Can this redistribution be done in such a way that equilibrium salaries are commensurate with the value of each scientist? Can we create an incentive-compatible academia, where it is in everyone's best interest to honestly report whose research they find valuable? Let me know what you think in comments.

## Wednesday, 25 January 2012

### Observe, don't just see, Dr Watson

The other day I was reading stories of Sherlock Holmes on the train. The following quite famous quote in which Holmes explains to Dr Watson his extraordinary inference and deduction skills caught my eyes:
Watson: And yet I believe that my eyes are as good as yours.
Holmes: Quite so. You see, but you do not observe...
In my interpretation this means that everybody is capable of seeing the same evidence, but what makes a good detective is the ability to focus on relevant aspects and putting these pieces together to reconstruct stories underlying their observations.

I started thinking about the relevance of this philosophy to data science and data-driven enterprises. Many companies, governments, etc now realise the value of data: the more data they have about their users/citizens, the better position they are in when making decisions - in theory. Take a look at for example the bazillion of startups based on the business model: "People like to engage in activity X. We will make their experience better by leveraging data from social networks". But can they really realise the potential in data? My impression is that many of them can't. These companies see; but they do not observe: they may be sitting on terabytes of data, but they're incapable of using it to tell the right story about their users.

For that, they need real data scientists, who are a scarce resource. Not every company can afford the luxury to have an experienced data science team playing with their data. Data scientists are private detectives of the modern world: they are presented with evidence (the data) and they are asked to uncover the hidden story that explains it all. In most cases, Sherlock Holmes had to find out who the murderer was and what their motivations were, data scientists (or rather, the algorithms they build) try to figure out what made you share a link on Facebook; whether your recent credit card transaction was made by you or it was fraud; or figure out from thousands of USB stick descriptions that "four gigabytes" and "4GB" mean the same thing.

For the average reader, the above examples would be somewhat less exciting than any of the murder cases Sir Arthur Conan Doyle's famous fictional consulting detective has ever worked on. But hey, maybe that's only because no-one with sufficient writing talent has picked this up yet: "Twelve Adventures of Hercule Bayes, the Consultant Data Scientist". I'd buy it.


text:
# What possible developments, theoretical or experimental, would convince you of the correctness or otherwise of each of the current approaches to quantum gravity?

There are several theoretical approach of vary degrees of success.

However, to test you want it to also describe theories we know, but the energy scales are so different how do we test them?

Eg. Particle colliders are still limited to the energy range of the standard model.

Cosmology seems perhaps more promising. So perhaps we need to look at quantum cosmology before we look for a theory of quantum gravity.

Gravitational wave detection from the early universe.

If we could make precise measurements of all over the universe, this could allow us to gain insight back to the big bang and early universe.

We could look for a dispersion relation of light in a vacuum.

For many people theoretical predictions will never been satisfying for promoting one theory over another, we really need some real experimental evidence.

There is not yet a single theory of quantum gravity. Throughout the years several approaches has been built in order to address the problem of quantizing gravity. In particular two current trends correspond to string theory and loop quantum gravity (LQG).  The developments, theoretical or experimental, that may prove or disprove the different theories about quantum gravity depend strongly on 'philosophical convictions'. That is,  there is no ultimate 'smoking gun' that will rule out all but one theory. The reason behind this is that most of the current approaches rely on different assumptions and, although they address the problem of quantizing gravity, their predictions may be different and may not overlap with one another. Thus, it is a matter of choice to determine which of such predictions will be considered as a proof of the validity of a theory of quantum gravity. Data coming form experiments like BICEP may  give some information about the so-called linearized quantum gravity. The data provided by the experiment serves as a motivation to seek for an underlying quantum theory of gravity. However, the results from this experiment may only offer information about a piece of the problem and may not serve as sufficient evidence to prove or disprove the current approaches. It has been suggested that an important step towards an experimental proof of a theory of quantum gravity is to observe a primordial black-hole in evaporation. This will at least give some data that can be compared to certain approaches to quantum gravity. One of the main requirements on a quantum theory of gravity is that it reduces to the standard general relativity in an appropriate limit. This at least serves a way of ruling out certain theories  since the description offered by general relativity agrees very well with the data obtained from experiments. In summary, the 'smoking gun' for a quantum theory of gravity is absent. There may be several future proposals that may suggest that one theory is better than the other, but this wont be sufficient to rule out other approaches. If all the approaches could provide the same observational predictions, then it will be an easy task to rule out some theories according to observations. In the current state of the field, this remains an open question without a general consensus.

The question was slightly reformulated: What developments would really bring you to drop a certain approach to quantum gravity.

As an example it was pointed out that any sign of Lorentz violating effect would drastically change the landscape of quantum gravity research but that so far it is surprising that no such effect has been found.

It was asked at which point should we drop SUSY as a credible explanation?

Any theory that would a posteriori reduce the number of (standard model) parameters would gain huge credibility. Likewise for explaining dark energy or dark matter.

Strings – smoking gun is extra dimensions/black holes

Is there any theory that can’t wiggle out of experimental evidence?

String theory – compare what it predicted in 1998 to what we see since then…

What would rule a theory out?  If it doesn’t have the right low energy limit.

Any experiment which show violation of quantum mechanics would kill most theories

What about violation of Lorentz invariance?

Is there a smoking gun signature for a bounce?

Good evidence if one can accommodate standard model type matter and predict a small positive cosmological constant

Possible metastability of the current vacuum – what happens after. Since it’s related to vacuum, also related to \lambda

In the case of string theory, it doesn’t seem that any experimental data could really convince us that the theory is true - or that it is not true. This is because string theory seems to accommodate every possible experimental result via the string landscape. It’s difficult to rule out a theory that is so flexible. Also, at this point string theory only describes perturbative quantum gravity; it can’t tell us anything about the black holes at the center of galaxies, for example.

Other approaches, such as loop quantum gravity, try to give a non-perturbative description of quantum gravity, which is harder to do but is more conceptually interesting. It could be that several of the existing approaches to quantum gravity will each contribute to a different part of the final theory, or a different regime of validity. Regarding loop quantum gravity, the question of whether we can obtain general relativity from it in the semiclassical limit still remains open. If this is achieved, it might convince some people that this theory is more worth working on.

[Bad Link: Plugin Not Found]


text:
# 100% Efficient Transfer of Energy

1. Sep 11, 2010

### grounded

Is there ANY way for ANY type of energy to be transformed from one form to another with no losses, a perfect 100% efficiency?

Thank you...

2. Sep 11, 2010

### ZapperZ

Staff Emeritus
Gravitational PE to KE. Drop a ball in low vacuum.

Zz.

3. Sep 11, 2010

### Staff: Mentor

Electrical to heat.

4. Sep 11, 2010

### grounded

I'm not arguing, just trying to understand if the 2nd law of thermodynamics applies only for heat to mechanical or if it applies to all forms off energy tranformation.

It seems logical to me that the amount of energy required to stop the ball from falling should always be more than just the energy needed to hold the weight of the ball.

My logic is that by suspending the ball, it will somewhat be compressed causing some amount of heat loss energy, although really small. ??

Also, I'm not to sure that PE turns into KE, it seems more like PE stops KE from increasing, also KE should increases with distance where PE should decrease with distance.

Any other examples?

5. Sep 11, 2010

### grounded

Wait, I'm wrong.
PE is not the energy required to hold the object from falling, its basically what KE can be at different distances. In that sense they are equal.

6. Sep 11, 2010

Its not possible to get 100% Efficient energy transfer, Only under Ideal condition it is possible, Ideal conditions are assumptions made by us.

In most of energy transfers taking place in the universe, losses take place generally due to friction, or the energy might be lost to atmosphere if its not a perfect closed system.

7. Sep 11, 2010

### ZapperZ

Staff Emeritus
They are? This doesn't make any sense. What if, instead of dropping it, I roll a wheel down an inclined plane. Do you think the KE of that wheel is the same as the dropped ball after it has gone through the same height? If it isn't, then how can PE be nothing more than "what KE can be at different distance" when these two examples produced DIFFERENT KE?

Zz.

8. Sep 11, 2010

### Dr Lots-o'watts

Energy is always transferred from one input type to many output types (the application + "losses") with perfect efficiency (law of conservation). The issue here is to single out one of the output types. It comes down to a matter of how close you get to what you actually want. It comes down to how perfect is heat isolation, how perfect is your crystal, how perfect is your vacuum etc. depending on what types of energy you work with.

9. Sep 11, 2010

### grounded

I must be missing something simple, I thought the PE of the ball at the top of the hill is equal to the KE of the ball at the bottom of the hill? (minus friction)

Isn't the potential amount of energy the ball can have equal to the the energy the ball can obtain by falling or rolling?

10. Sep 11, 2010

### ZapperZ

Staff Emeritus
Not if the ball is rolling! The original PE has been transferred into the KE of the ball and the rotational energy.

Still, this is besides the point. By definition, the potential energy is never defined as what it's KE can be. The potential energy for an electron in a stable state in an atom has nothing to do with what its KE can be when it is liberated from that atom.

Zz.

11. Sep 11, 2010

### Staff: Mentor

Look at the definition of change in entropy, as it appears in most textbooks:

$$dS = \frac{dQ}{T}$$

If dQ = 0, then dS = 0.

12. Sep 11, 2010

### Bob S

If you look at the transfer of energy from one nuclear isomer to another identical nucleus by way of an emitted photon, the energy has to be conserved. Look in particular at Fe-57 isomer used in Mossbauer experiments. The emitted photon energy spread is
~5 *10-9 eV compared to the emitted photon energy of 14.4 *103 eV, so the energy transfer has to be 100% to within 1 part in 1012 to absorbed by another Fe-57 nucleus. The major energy loss is due to nuclear recoil because of the finite nuclear mass. This can be minimized by 'locking" the Mossbauer atoms in a crystal lattice. (Sorry, this is not Classical Physics.)

Bob S


text:
# Math Help - graphing area

2. ## Re: graphing area

Hey vegasgunner.

Hint: Calculate the area of the triangle and use symmetry. Alternatively you can setup the appropriate integral and obtain the same answer.


text:
# Parameter identification for first order ODE

I have two arrays $f(z)$ and $z$ both indexed by k and I want to solve $\frac{df}{dz}=\mu(1-f)$ to find $\mu(z)$

What would be the best numerical method to solve this equation?

• this is a ODE parameter identification problem: given the two data arrays $(z_k,f_k)$ (indexed by $k$), find the coefficient $\mu$ that best fits the data, where best is defined in a statistically precise way. My first approach to this would be to write it as an ODE-constrained optimization problem with a properly regularized cost function, then solve its discrete equivalent using a numerical method. The slides I linked to should get you started if you decide to go this route. – GoHokies Aug 29 '17 at 14:24
• alternatively, MATLAB can help with this as well – GoHokies Aug 29 '17 at 14:24
• or, since you have mentioned FORTRAN, there is the PDEFIT package that seems to do the trick: paper and software – GoHokies Aug 29 '17 at 14:36
• I'm not sure about that recommendation in the slides. For nonlinear ODEs, solving parameter estimation problems with Newton's method or other local optimizers is usually a bad idea and will not find the local minima (though this case may be simple enough). Here I wrote up a general overview of how the algorithms tend to work though, and you can stick any optimizer and diffeq solver of your choice into that general framework. Of course, Julia has this all setup, but MATLAB has AMIGO2 as well – Chris Rackauckas Aug 29 '17 at 18:54
• @ChrisRackauckas usually a bad idea - that's a bold statement. Newton will usually converge to a local minimum. The "quality" of that minimum will depend on your prior, on the data quality ($p(\rm{data}|\mu)$) and the exact choice of local optimizer, among other factors. – GoHokies Aug 29 '17 at 20:00

What you have is a linear ODE-constrained parameter identification problem: given the two data arrays $(z_k,f_k)$ (indexed by $k$) and the ODE model (the constraint), find the parameter $\mu$ that best fits the data, where the notion of best is defined in a statistically precise way. A very simple example would be a least-squares type regularized functional

$$\arg \min_\mu \sum_k \left( f_k - f(z_k) \right)^2 + \| \mu \|^2 \; {\rm subject \; to\;} f_z = \mu (1-f).$$

One approach to this would be to write it as an ODE-constrained optimization problem with a properly regularized cost function, then solve its discrete equivalent using a numerical method and off-the shelf optimizer, as exemplified in these slides (1). Note that the optimization algorithm used there is a local Newton-type procedure. This means, for most practical problems, that you will need a reasonably good initial guess for the Newton procedure to converge to a "good-quality" minimum point, i.e., to a $\bar\mu$ not "too far" (in some norm) from the unknown "true" parameter $\mu^{\rm true}$.

If local minima are a major concern, and you do not have a good "guess" at $\mu^{\rm true}$ to start your Newton routine from, then global optimization is an alternative (caveat emptor: these procedures are more computationally expensive than a local routine). Or you can bring out the best of the two approaches into a combined global/local optimization strategy.

If you don't have the time to implement all this, then there are other options to consider:

1. The MATLAB system identification toolbox can help (Google idgrey, greyest). The workflow for your problem is described in detail here.
2. There are some FORTRAN libraries that do this. For example PDEFIT and/or EASY-FIT. I'd suggest you to read the manual to get acquainted with the features (and limitations) of each.
3. Other programming languages appear to have at least some support for parameter identification: R, python with numpy/scipy.


text:
## Distributed robust $$H_\infty$$ consensus control in directed networks of agents with time-delay.(English)Zbl 1140.93355

Summary: This paper investigates consensus problems for directed networks of agents with external disturbances and model uncertainty on fixed and switching topologies. Both networks with and without time-delay are taken into consideration. In doing the analysis, we first perform a model transformation and turn the original system into a reduced-order system. Based on this reduced-order system, we then present conditions under which all agents reach consensus with the desired $$H_\infty$$ performance. Finally, simulation results are provided to demonstrate the effectiveness of our theoretical results.

### MSC:

 93B36 $$H^\infty$$-control 93A15 Large-scale systems 93C41 Control/observation systems with incomplete information 93B11 System structure simplification
Full Text:

### References:

 [1] Tsitsiklis, J.N.; Bertsekas, D.P.; Athans, M., Distributed asynchronous deterministic and stochastic gradient optimization algorithms, IEEE trans. automat. control, 31, 9, 803-812, (1986) · Zbl 0602.90120 [2] Bertsekas, D.P.; Tsitsiklis, J., Parallel and distributed computation, (1989), Prentice-Hall Upper Saddle River, NJ · Zbl 0743.65107 [3] Lynch, N.A., Distributed algorithms, (1997), Morgan Kaufmann San Mateo, CA [4] Vicsek, T.; Cziroók, A.; Ben-Jacob, E.; Cohen, O.; Shochet, I., Novel type of phase transition in a system of self-driven particles, Phys. rev. lett., 75, 6, 1226-1229, (1995) [5] Jadbabaie, A.; Lin, J.; Morse, A.S., Coordination of groups of mobile autonomous agents using nearest neighbor rules, IEEE trans. automat. control, 48, 6, 988-1001, (2003) · Zbl 1364.93514 [6] Moreau, L., Stability of multi-agent systems with time-dependent communication links, IEEE trans. automat. control, 50, 2, 169-182, (2005) · Zbl 1365.93268 [7] Olfati-Saber, R.; Murray, R.M., Consensus problems in networks of agents with switching topology and time-delays, IEEE trans. automat. control, 49, 9, 1520-1533, (2004) · Zbl 1365.93301 [8] Olfati-Saber, R.; Fax, J.A.; Murray, R.M., Consensus and cooperation in networked multi-agent systems, Proc. IEEE, 95, 1, 215-233, (2007) · Zbl 1376.68138 [9] Fax, J.A.; Murray, R.M., Information flow and cooperative control of vehicle formations, IEEE trans. automat. control, 49, 9, 1465-1476, (2004) · Zbl 1365.90056 [10] Lafferriere, G.; Williams, A; Caughman, J.; Veerman, J.J.P., Decentralized control of vehicle formations, System control lett., 54, 9, 899-910, (2005) · Zbl 1129.93303 [11] D.B. Kingston, R.W. Beard, Discrete-time average-consensus under switching network topologies, in: Proceedings of IEEE Conference on Decision and Control, 2005, pp. 2016-2021 [12] Ren, W.; Beard, R.W., Consensus seeking in multi-agent systems under dynamically changing interaction topologies, IEEE trans. automat. control, 50, 5, 655-661, (2005) · Zbl 1365.93302 [13] W. Ren, R.W. Beard, E.M. Atkins, A survey of consensus problem in multiagent coordination, in: Proceedings of the American Control Conference, 2005, pp. 1859-1864 [14] Olfati-Saber, R., Flocking for multi-agent dynamic systems: algorithms and theory, IEEE trans. automat. control, 51, 3, 410-420, (2006) · Zbl 1366.93391 [15] H. Tanner, A. Jadbabaie, G. Pappas, Stable flocking of mobile agents, part II: Dynamic topology, in: Proceedings of IEEE Conference on Decision and Control, 2003, pp. 2016-2021 [16] Moshtagh, N.; Jadbabaie, A., Distributed geodesic control laws for flocking of nonholonomic agents, IEEE trans. automat. control., 52, 4, 681-686, (2007) · Zbl 1366.93390 [17] Dimarogonas, D.V.; Kyriakopoulos, Kostas J., On the rendezvous problem for multiple nonholonomic agents, IEEE trans. automat. control, 52, 5, 916-922, (2007) · Zbl 1366.93401 [18] D.V. Dimarogonas, Kostas J. Kyriakopoulos, Connectivity preserving distributed swarm aggregation for multiple kinematic agents, in: Proceedings of IEEE Conference on Decision and Control, 2007, pp. 2913-2918 [19] M. Ji, M. Egerstedt, Distributed formation control while preserving connectedness, in: Proceedings of IEEE Conference on Decision and Control, 2006, pp. 5962-5967 [20] Hong, Y.; Gao, L.; Cheng, D.; Hu, J., Tracking control for multi-agent consensus with an active leader and variable topology, Automatica, 42, 7, 1177-1182, (2006) · Zbl 1117.93300 [21] Ögren, P.; Fiorelli, E.; Leonard, N.E., Cooperative control of mobile sensor networks: adaptive gradient climbing in a distributed environment, IEEE trans. automat. control, 49, 8, 1292-1302, (2004) · Zbl 1365.93243 [22] Xiao, L.; Boyd, S., Fast linear iterations for distributed averaging, System control lett., 53, 1, 65-78, (2004) · Zbl 1157.90347 [23] Xiao, L.; Boyd, S.; Kimb, S.J., Distributed average consensus with least-Mean-square deviation, J. parallel distr. comput., 67, 1, 33-46, (2007) · Zbl 1109.68019 [24] P. Bliman, G. Ferrari-Trecate, Average consensus problems in networks of agents with delayed communications, in: Proceedings of IEEE Conference on Decision and Control, 2005, pp. 7066-7071 [25] L. Fang, P.J. Antsaklis, Information consensus of asynchronous discrete-time multi-agent systems, in: Proceedings of the American Control Conference, 2005, pp. 1883-1888 [26] M. Cao, A.S. Morse, B.D.O. Anderson, Agreeing asynchronously: Announcement of results, in: Proceedings of IEEE Conference on Decision and Control, 2006, pp. 4301-4306 [27] M. Cao, A.S. Morse, B.D.O. Anderson, Reaching an agreement using delayed information, in: Proceedings of IEEE Conference on Decision and Control, 2006, pp. 3375-3380 [28] Kim, Y.; Mesbahi, M., On maximizing the second smallest eigenvalue of a state-dependent graph Laplacian, IEEE trans. automat. control, 51, 1, 116-120, (2007) · Zbl 1366.05069 [29] P. Lin, Y. Jia, J. Du, S. Yuan, Distributed consensus control for networks of second-order agents with fixed topology and time-delay, in: Proceedings of the Chinese Control Conference, 2007, pp. 2-577-2-581 [30] Lin, P.; Jia, Y., Average-consensus in networks of multi-agents with both switching topology and coupling time-delay, Physica A, 387, 1, 303-313, (2008) [31] P. Lin, Y. Jia, J. Du, S. Yuan, Distributed control of multi-agent systems with second-order agent dynamics and delay-dependent communications, Asian J. Control 10 (2) (2008) (in press) [32] Godsil, C.; Royle, G., Algebraic graph theory, (2001), Springer-Verlag New York · Zbl 0968.05002 [33] Horn, R.A.; Johnson, C.R., Matrix analysis, (1987), Cambridge Univ. Press Cambridge [34] Dullerud, G.E.; Paganini, F., () [35] K.C. Goh, L. Turan, M.G. Safonov, G.P. Papavassilopoulos, J.H. Ly, Baffine matrix inequality properties and computational methods, in: Proceedings of the American Control Conference, 1994, pp. 850-855 [36] S. Ibaraki, M. Tomizuka, Rank minimization approach for solving BMI problems with random search, in: Proceedings of the American Control Conference, 2001, pp. 1870-1875 [37] Kanev, S.; Scherer, C.; Verhaegen, M.; De Schutter, B., Robust output-feedback controller design via local BMI optimization, Automatica, 40, 7, 1115-1127, (2004) · Zbl 1051.93042 [38] Sastry, S., Nonlinear systems: analysis, stability, and control, (1999), Springer [39] Hale, J., Theory of functional differential equations, (1977), Springer-Verlag New York
This reference list is based on information provided by the publisher or from digital mathematics libraries. Its items are heuristically matched to zbMATH identifiers and may contain data conversion errors. It attempts to reflect the references listed in the original paper as accurately as possible without claiming the completeness or perfect precision of the matching.