text:
Parameters: Simulation

Command ElementModifies the simulation parameter element.

Format

<Param_Simulation
[ constr_tol          = "real" ]
[ implicit_diff_tol   = "real" ]
[ nr_converge_ratio   = "real" ]
[ linsolver           = "HARWELL" ]
[ acusolve_cosim      = { "TRUE" | "FALSE" } ]
[ contact_iteration   = { "ACTIVE" | "INACTIVE" } ]
/>

Attributes

constr_tol
Modifies the accuracy to which the system configuration and motion constraints are to be satisfied at each step. This is used by kinematics and the transient solver in ODE formulation (ABAM, VSTIFF, or MSTIFF), but not the DAE form, which has its own tolerance for this (Param_Transient::dae_constr_tol). This should be a small positive number.
The default value for constr_tol = 1.0E-10.
implicit_diff_tol
Modifies the accuracy to which implicit differential equations, such as Control_Diff equations with the is_implicit = "TRUE", are to be satisfied. The default value for implicit_diff_tol = 1.0E-5.
nr_converg_ratio
Modifies a measure of the rate of convergence in the Newton-Raphson method for ODE solvers. If the maximum entry in the constraint vector is larger than nr_converg_ratio times the maximum entry from the previous iteration, the Newton-Raphson iterations are converging slowly and the generalized coordinates will be re-partitioned in the next integration step to select a new set of independent coordinates. This attribute is applicable only if an ODE integrator (ABAM, VSTIFF, or MSTIFF) is selected.
The default value is 0.09.
linsolver
The type of linear solver used, which currently is set to "HARWELL" for all analyses.
acusolve_cosim
A Boolean flag that determines whether the simulation will be coupled with AcuSolve or not.
"TRUE" means that the MBD model is to be coupled with the CFD model in a co-simulation between MotionSolve and AcuSolve.
"FALSE" means that the MBD model is not to be coupled with the CFD model.
The default value is "FALSE".
contact_iteration
Specifies whether the contact residual is being updated during the corrector step of the implicit solver and an analytical Jacobian is provided or not (see Comment 2).

The value "ACTIVE" indicates that the residual is being updated. The value "INACTIVE" indicates that the residual is not being updated.

The attribute contact_iteration is optional. When not specified, it defaults to "INACTIVE".

Example

The example below shows the default settings for the Param_Simulation command element.

<Param_Simulation
constr_tol          = "1E-10"
implicit_diff_tol   = "1E-5"
nr_converg_ratio    = "0.09"
linsolver           = "HARWELL"
/>

1. Linear Solvers: At the core of most analyses in MotionSolve is solving a set of linear equations (A x = b). For example, the Newton-Raphson method solves a set of linear equations as part of the iteration process to find the solution to a set of non-linear algebraic equations. A brief explanation of how the linear solver works follows:
• The Harwell Linear Solver is a tool for solving linear algebraic equations. It is especially suited for linear systems characterized by non-singular, unsymmetrical and sparse coefficient matrices A that have a fixed, non-zero entry pattern. The latter implies that while matrix entries are allowed to change with time, the pattern of non-zeroes is not. The matrix entries must be real.
This methodology is suitable for solving small to medium sets of equations (<10,000 equations) and is therefore quite suitable for multi-body systems simulation. The Harwell software solves the equations in three major steps:
1. Symbolic LU factorization.
2. Numeric LU factorization.
3. Forward-Backward Substitution.
Symbolic LU factorization
Given a pattern of non-zero entries in A and their representative values, this function computes the symbolic lower- and upper-triangular (LU) factors of A. A partial pivoting scheme is used. It tries to maximize the stability of the LU factors while still maintaining sparsity of the factors. This operation is typically done once or only a few times during the simulation.
Numeric LU factorization
Given the current values of the non-zero entries in A and the symbolic LU factors, this utility returns the numeric LU factors of A. The symbolic LU factorization must therefore precede the numeric LU factorization. This operation is done whenever a new Jacobian is needed.
Forward-Backward Substitution
Given the numeric LU factors of a sparse matrix A and an appropriately sized right-hand-side vector b, this utility returns the solution x for the linear problem. The Numeric LU factorization must precede the forward-backward substitution operation. This operation is performed at each iteration.
2. The feature contact_iteration can increase the simulation robustness for contacting bodies, especially if there are contacts connecting multiple bodies or if high oscillations are observed. In certain circumstances, it can also decrease the computational performance. When using active contact iteration, contacts can react more sensitively when switching from an active to an inactive state, especially if the contacting bodies carry a high velocity difference. In those scenarios, using an additional zero-crossing sensor is recommended.


text:
by using Biot Savart Law, the calculation is very terrible. since it involve the cubic inverse distance. However, using vector potential is a way to make an analytic solution.  first in cylindrical coordinate, the current density is:

$\displaystyle \vec{J} = J_\phi \hat{\phi} = I \delta(z) \delta(\rho-a) \hat{\phi}$

the vector potential, without loss of generality, we set the field point on x-z plane and the vector potential is pointing $\hat{\phi}$ direction.

Oh, shit, i am tired on typing equation by Latex. anyway,

the finial result can be obtained in J.D. Jackson’s book. and it take the form of Elliptic Integral.

$\displaystyle A_\phi = \frac{\mu_0 I}{2\pi } \sqrt{\frac{a}{\rho}}\left( \frac{k^2-2}{2k} K(k^2) + \frac{1}{k} E(k^2) \right)$

where K(m) is the complete elliptic integral of first kind, and E(m) is the complete elliptic integral of second kind.

$\displaystyle K(m) = \int_0^{\pi/2} \frac{1}{\sqrt{1-m \sin^2(\theta)}}d\theta$

$\displaystyle E(m) = \int_0^{\pi/2} \sqrt{1-m \sin^2(\theta)} d\theta$

and

$\displaystyle k^2 = \frac{ 4 a \rho }{ (a+\rho)^2+z^2 }$

and using the derivative properties of elliptic integral, we can have the magnetic field in analytic form. The coil radius is $a$, at x-y plan. The observation point is $(\rho, z)$

$\displaystyle B_\rho(\rho,z) = \frac{\mu_0 I}{2\pi}\frac{z}{\rho \sqrt{(a+\rho)^2+z^2}} \left(K(k^2) - \frac{a^2+\rho^2 + z^2}{(a-\rho)^2 + z^2} E(k^2) \right)$

$\displaystyle B_z(\rho,z) = \frac{\mu_0 I}{2\pi}\frac{1}{ \sqrt{(a+\rho)^2+z^2}} \left(K(k^2) + \frac{a^2-\rho^2 - z^2}{(a-\rho)^2 + z^2} E(k^2) \right)$

Here is the PDF for the detail calculation. Field of Single Coil. Note that, in the pdf, the B-field is smaller by factor of 2. The denominator should be $2\pi$. As we check the magnetic field on the z-axis, the above formula will reduce to

$\displaystyle B_z(0,z) = \frac{\mu_0 I}{2}\frac{a^2}{ (a^2+z^2)^{3/2}}$

where $K(0) = E(0) = \pi/2$

format: plain


text:
# EE 381 Assignment 06 – Markov Chains solution

\$25.00

Category:

## Description

0. Introduction and Background Material
0.1. A simple example of a Markov chain
A machine can be in two possible states A and B. The state of the machine at time ( 1) k + is
determined by the state of the machine at the previous time step k :
• If the state at k is A, then the state at ( 1) k + will be A with probability 11 p
• If the state at k is A, then the state at ( 1) k + will be B with probability 12 p
• If the state at k is B, then the state at ( 1) k + will be A with probability 21 p
• If the state at k is B, then the state at ( 1) k + will be B with probability 22 p
These state transitions are shown schematically in Figure 1. Also, Figure 2 shows the same
information on state transition probabilities in a more compact form.

Figure 0.1 Figure 0.2
Clearly: 11 12 p p + =1 and 21 22 p p + =1.
A p11 p12
TIME k TIME (k+1)
B
A
B
p22
p21
A
B
p12
p21
p11
p22
EE 381 Project : Markov Chains Dr. Chassiakos – Page 2
For the current problem the transition probabilities have the following values:
11 12 21 22 pppp = 0.8 ; 0.2 ; 0.5 ; 0.5 = = =
The initial state is randomly distributed between state “A” and state “B”.
The probability of the initial state being “A” is given as: Prob (Initial State = A) = 0.4
The probability of the initial state being “B” is given as: Prob (Initial State = B) = 0.6
The graphs below show several simulation results for this two-state Markov chain. The chain is
run for n =15 time steps.
Figure 0.3 shows a single simulation run of the chain. The initial state of the chain is “A”. The
graph shows the evolution of the chain for the next 15 steps. In this particular simulation the
chain goes through the following sequence of states, after the initial state “A”:
[A BBABBAABBBAAABA]
Figure 0.4 shows another single simulation run of the chain. In this case the initial state of the
chain is “B”. The graph shows the evolution of the chain for the next 15 steps. In this particular
simulation the chain goes through the following sequence of states, after the initial state “B”:
[B AAABAAAAAAAAABA]
Figure 0.3 Figure 0.4
EE 381 Project : Markov Chains Dr. Chassiakos – Page 3
Figure 0.5 shows the combined results from N=10000 simulation runs of the chain. It is seen that
the chain starts at state “A” 40% of the time, and at state “B” 60% of the time. The graph shows
the evolution of the chain for the next 15 steps.
Figure 0.6 shows the results of the state evolution using the “State Transition Matrix” approach.
The probabilities are CALCULATED based on the transition matrix. The probabilities they ARE
NOT derived as the result of N=10000 simulation runs, as was done in Figure 0.5. It is noted that
the results of Figures 0.5 and 0.6 are very similar, almost identical, confirming the fact that the
state transition matrix approach can be used instead of a step-by-step simulation.
Figure 0.5 Figure 0.6
EE 381 Project : Markov Chains Dr. Chassiakos – Page 4
1. A three-state Markov Chain
Follow the previous code as a guide, and simulate the Markov chain defined in problem 11.1 of
the AMS text by Grinstead & Snell, p. 406 (reference [1]).
The transition probabilities are given by the following matrix:
111
333
1 1 0
2 2
111
442
P
 
 
 
  =  
 
     
1. The initial probability distribution of the state is given by: [ ] 111
424
RNS   =    
2. Run each experiment for n =15 time steps. In order to obtain meaningful statistical data
perform a total of N =10,000 experiments.
3. After your experimental simulations are completed, use the State Transition Matrix approach
and compare the results.
SUBMIT a report with the results and your code. You must follow the guidelines given in the
syllabus regarding the structure of the report. Points will be taken off, if you do not follow the
guidelines. The report should contain:
• One plot showing one single-simulation run for n =15 steps, similar to Figure 0.3, but for a
three-state chain. For an example see the figure below. Note that this figure does not have
labels, but your plot should have the appropriate labels and title.
• A plot with the experimental probabilities for the states R, N, S, similar to Figure 0.5
• A plot with the probabilities obtained through the state transition matrix approach for the
three states R, N, S, similar to Figure 0.6
• The code in an Appendix
• Make sure that all the plots are properly labeled
EE 381 Project : Markov Chains Dr. Chassiakos – Page 5
This problem presents an introduction to the algorithm used for ranking web pages for searching
purposes. The algorithm was developed by Google’s founders Page and Brin with Motwani and
Winograd, and was first published in 1998 (see reference [2]). The PageRank algorithm allowed
Google to rise to the top of all web search engines within a matter of months after its
implementation, outperforming the established search engines of the time such as AltaVista and
Yahoo. The algorithm is based on the theory of Markov chains, utilizing the information on the
The current problem uses a simplified web in order to show how the algorithm works (see
reference [3]). The simplified web consists of 5 pages only, and it is shown schematically in
Figure 2.1.
Figure 2.1: A five-page web
A Markov chain model of an impartial surfer visiting web pages is created as following:
• At time k the surfer is on page X , where X ABCDE ∈{ ,,,, } .
• At time ( 1) k + the surfer will move randomly to another page Y , which must be linked to X .
• The state k S of the Markov chain at time k is defined as the page which the surfer is visiting
at time k : S ABCDE k ∈{ ,,,, } .
• To create the Markov chain model, the algorithm assumes that all pages that can be reached
from X have the same probability to be visited by the surfer. Thus:
( 1 )
1 , if can be reached from (Total number of links leaving page X)
Prob |
0 , if cannot be reached from
k k
Y X
S YS X
Y X
+

 = = = 


• For example: 1
1 Prob ( | ) 3 k k S AS C + = = = ; Prob ( | ) 0 k k 1 S ES D + = = = ; etc.
A
B
D
C
E
EE 381 Project : Markov Chains Dr. Chassiakos – Page 6
• Based on these probabilities the State Transition Matrix ( P ) of the Markov chain is
constructed. The left eigenvector w of P corresponding to the eigenvalue λ =1 (which is a
fixed vector of the Markov chain) is computed from: wP w=
• The PageRank algorithm uses the vector w to rank the pages of the five-page web in terms of
visiting importance.
To complete this problem follow the steps below:
1. Create the 5 5 × State Transition Matrix P with the values of the transition probabilities
2. Assume that initially all pages are equally likely to be visited, i.e. use the initial state
probability distribution vector: 1 [ ] 11111
55555
v ABCDE   = =    
3. Calculate the probability vectors for n =1,2, 20  steps, using the State Transition
matrix only. Note that you ARE NOT asked to run the complete simulation of the chain. You
are only asked to use the State Transition matrix approach to calculate the probability
vectors for n steps.
4. Rank the pages {ABCDE ,,,, } in order of importance based on the results from the previous
step
5. Create a plot of the calculated state probabilities for each of the five states {ABCDE ,,,, } vs.
the number of steps for n =1,2, 20  . This should be similar to the plot in Figure 0.6
6. Assume that the surfer always starts at his/her home page, which for the purposes of this
problem is page E . Repeat steps (2)-(5) with the new state probability distribution vector:
v ABCDE 2 = [ ] = [00001]
SUBMIT a report with the results and your code. You must follow the guidelines given in the
syllabus regarding the structure of the report. Points will be taken off, if you do not follow the
guidelines. The report should contain:
• The 5 5 × State Transition Matrix P
• For the initial state probability distribution vector 1 v :
 Submit the page ranking calculated in Step 4. Use Table 1 for your answer. Note: You will
need to replicate the table in your Word file, in order to provide the answer in your report.
Points will be taken off if you do not use the table.
 Submit the plot created in Step 5.
• For the initial state probability distribution vector 2 v
 Submit the page ranking calculated in Step 4. Use Table 1 for your answer. Note: You will
need to replicate the table in your Word file, in order to provide the answer in your report.
Points will be taken off if you do not use the table.
 Submit the plot created in Step 5.
• The code in an Appendix.
• Make sure the plots are properly labeled

EE 381 Project : Markov Chains Dr. Chassiakos – Page 7
Initial probability vector: 1 v Initial probability vector: 2 v
Rank Page Probability Rank Page Probability
1 1
2 2
3 3
4 4
5 5
Table 1.
References
[1] “Introduction to Probability”, C.M. Grinstead and J.L Snell, American Mathematical
Society, electronic version, 2003, 2nd edition
[2] L. Page, S. Brin, R. Motwani, and T. Winograd. “The PageRank Citation Ranking:
Bringing Order to the Web”, Technical Report, Stanford University, 1998.
[3] “Mathematics & Technology”, C. Rousseau & Y. Saint-Aubin, Springer, 2008.


text:
Question

# Three containers have their volumes in the ratio 3:4:5. They are full of milk and water. The mixtures contain milk and water in the ratios of 4:1, 3:1 and 5:2, respectively. The ratio of milk and water in the fourth container which contains all these mixtures is [a] 4:1[b] 151:48[c] 157:53[d] 5:2

Hint: Assume that the volumes of the containers are 3x,4x and 5x. Find the content of milk and water in each of these solutions. Hence find the content of milk and water in the final solution. Hence find the ratio of the content of milk to the content of water.

Complete step-by-step solution -

Let the volume of the first container be 3x, the volume of the second container be 4x, and that of the third container be 5x.
Amount of milk in the first container $=\dfrac{4}{1+4}3x=\dfrac{12x}{5}$
Amount of milk in the second container $=\dfrac{3}{3+1}4x=\dfrac{12x}{4}=3x$
Amount of milk in the third container $=\dfrac{5}{5+2}5x=\dfrac{25x}{7}$
Amount of water in the first container $=\dfrac{1}{4+1}3x=\dfrac{3x}{5}$
Amount of water in the second container $=\dfrac{1}{3+1}4x=x$
Amount of water in the third container $=\dfrac{2}{2+5}5x=\dfrac{10x}{7}$
Hence the total amount of milk in the final container $=\dfrac{12x}{5}+3x+\dfrac{25x}{7}=\dfrac{84x+105x+125x}{35}=\dfrac{314x}{35}$
The total amount of water in the final container $=\dfrac{3x}{5}+x+\dfrac{10x}{7}=\dfrac{21x+35x+50x}{35}=\dfrac{106x}{35}$
Hence the ratio of the amount of milk to the amount of water in the fourth container $=\dfrac{314x}{35}:\dfrac{106x}{35}=314:106=157:53$
Hence option [c] is correct.

Note: Alternatively, you can find the ratio of milk to the total volume and hence find the ratio of milk to water.
Let the total volume of the final solution be (3+4+5)=12.
So we have the amount from the first solution = 3
Hence milk from the first solution $=\dfrac{12}{5}$
Similarly milk from the second solution $=3$ and the milk from the third solution $=\dfrac{25}{7}$
Hence total milk $=\dfrac{12}{5}+3+\dfrac{25}{7}=\dfrac{314}{35}$
Hence the ratio of milk to the total volume $=\dfrac{314}{35}:12=\dfrac{157}{35}:6=157:210$
Hence the ratio of milk to water $=157:\left( 210-157 \right)=157:53$


text:
# Reply to 'Flawed assumptions compromise water yield assessment'

## Article metrics

• 1094 Accesses

• 1 Citations

## Introduction

We appreciate the comments on the relative contributions of climate and watershed condition to water yield by Gudmundsson et al.1 (hereafter G17) on our study2 (hereafter Z15), both of which investigated the global pattern of the effect of climate and land cover on water yield. We acknowledge their adoption of a more commonly used definition on the relative contribution. In this correspondence, we confirm the validity of the major conclusions of Z15 and compare the definitions and calculations of the relative contributions between Z15 and G17. We also identify an important research gap on the most appropriate estimation of the relative contributions of climate and land cover to the annual water yield using these global theoretic frameworks (e.g., Fuh3 and Budyko), which require further attention.

Z15 first validated the broad application of the Fuh framework (Eqs. 1 and 2) using published watershed data across the globe. Then, Z15 used the theoretical reasoning of the Fuh framework to identify two critical values (ψ = P/PET = 1 and m = 2) and their corresponding hydrological sensitivities and concluded that changes in land cover in watersheds with ψ < 1 or m < 2 can lead to greater hydrological alterations. In addition, Z15 applied those results to estimate the relative contributions of ψ and m to R/P, which pertain to an application aspect of the key conclusions.

G17 questioned the calculation of the relative contributions of Z15. It is critical to stress that such questioning was only on the application aspect of Z15 rather than on their central conclusions. However, the title of G17 has led to significant confusion and a misunderstanding of Z15 as well as of global theoretical frameworks (Fuh, Budyko, etc.) in scientific communities. Thus, we further confirm that the central conclusions of Z15 were not questioned by G17.

G17 also raised a question regarding the physical inconsistency between the results shown in Fig. 7a of Z15 and the central statements presented in the abstract of Z15. The results presented in Fig. 7a of Z15 are about the relative contributions (Eq. (3)), whereas the central statements in the abstract of Z15 refer to sensitivity (Eqs. 10 and 11). Here, we acknowledge that “more responsive” in the abstract of Z15 should be “more sensitive,” and this imprecise wording may lead to G17’s misunderstanding of Fig. 7a of Z15 and the statement in the abstract of Z15.

$$\frac{R}{P} = Z\left( {\psi ,m} \right) = \left( {1 + \psi ^{ - m}} \right)^{1/m} - \psi ^{ - 1}{,}$$
(1)

or

$$\frac{R}{P} = Z\left( {\phi ,m} \right) = \left( {1 + \phi ^m} \right)^{1/m} - \phi{,}$$
(2)

where ψ is defined as the wetness index ($$\psi = P/{\mathrm{PET}}$$ in Eq. (1) and ϕ is the dryness index ($$\phi = {\mathrm{PET}}/P$$ in Eq. (2)).

The central debate between Z15 and G17 lies in the relative contributions of ψ and m to R/P. The relative contributions can be defined in various ways depending on the research objectives. Z15 defined them as the relative contributions of ψ and m to the sensitivities of R/P (Eq. (3)), aiming to compare the relative magnitudes of the sensitivity of ψ and m to the sensitivities of R/P. By contrast, G17 adapted a more commonly used definition, as shown in Eq. (5), which can be interpreted as the changed amounts of R/P due to the changes in ψ and m. Clearly, two definitions reflect different physical variables of interest, and consequently they are not comparable.

$$C_{m,\psi } = \frac{{100}}{{1 + \left| \alpha \right|}},\,{\mathrm{where}} \, \alpha = \left( {\frac{{\partial Z}}{{\partial \psi }}} \right)/\left( {\frac{{\partial Z}}{{\partial m}}} \right){,}$$
(3)
$$C^ \ast _{m,\phi } = \frac{{100}}{{1 + \left| {\alpha ^ \ast } \right|}},\,{\mathrm{where}} \, \alpha ^ \ast = \left( {\frac{{\partial Z^ \ast }}{{\partial \phi }}} \right)/\left( {\frac{{\partial Z^ \ast }}{{\partial m}}} \right){,}$$
(4)
$$C_{m,corr} = \frac{{100}}{{1 + \left| \beta \right|}},\,{\mathrm{where}} \, \beta = \left( {\frac{{\partial Z}}{{\partial \psi }}\Delta \psi } \right)/\left( {\frac{{\partial Z}}{{\partial m}}\Delta m} \right){.}$$
(5)

Using the definition of Z15, G17 further challenged the physical consistency of the approach suggested by Z15 due to the rearrangement of the Fuh framework from ψ (P/PET) to ϕ (PET/P) (Eqs. (3) and (4)). As a result, G17 found two different patterns (Fig. 1a, c in G17). We believe that G17 made a simple mistake by assuming Cm (Eq. (3)) and C*m (Eq. (4)) to be the same. Although they all refer to the relative contributions of m, the former is relative to the climate with regard to wetness (P/PET), whereas the latter is relative to the climate with regard to dryness (PET/P). Therefore, it is no surprise that they have different trends.

To further demonstrate the physical consistency of Z15, Eqs. (3) and (4) can be rewritten as Eqs. (6) and (7), respectively, when ϕ is rearranged to ψ (or ψ to ϕ; see Eqs. (8)–(10) for the derivation in detail). As shown in Figs. 1 and 2, Cm,ψ and C, are identical in disregard to the rearrangement from ψ (P/PET) to ϕ (PET/P). The same phenomenon is also found for C* and C*. Therefore, Figs. 1 and 2 clearly prove that the approaches in Z15 are physically consistent.

$$C_{m,\phi } = \frac{{100}}{{1 + \left| {\frac{{\partial Z}}{{\partial \phi }} \times \frac{{\mathrm{d}\phi }}{{\mathrm{d}\psi }}{\mathrm{/}}\frac{{\partial Z}}{{\partial m}}} \right|}} = \frac{{100}}{{1 + \left| { - \frac{1}{{\psi ^2}} \times \frac{{\partial Z^ \ast }}{{\partial \phi }}{\mathrm{/}}\frac{{\partial Z^ \ast }}{{\partial m}}} \right|}} = \frac{{100}}{{1 + \left| {\phi ^2\alpha ^ \ast } \right|}}{,}$$
(6)
$$C^ \ast _{m,\psi } = \frac{{100}}{{1 + \left| {\frac{{\partial Z^ \ast }}{{\partial \psi }} \times \frac{{\mathrm{d}\psi }}{{\mathrm{d}\phi }}{\mathrm{/}}\frac{{\partial Z^ \ast }}{{\partial m}}} \right|}} = \frac{{100}}{{1 + \left| { - \frac{1}{{\phi ^2}} \times \frac{{\partial Z}}{{\partial \psi }}{\mathrm{/}}\frac{{\partial Z}}{{\partial m}}} \right|}} = \frac{{100}}{{1 + \left| {\psi ^2\alpha } \right|}}{.}$$
(7)

## Methods

### Sensitivities of R/P to ψ (Sψ) and m (Sm)

The sensitivities of R/P to ψ (Sψ) and m (Sm) can be derived from Eq. (1) as:

$$S_\psi = \frac{{\partial Z}}{{\partial \psi }} = \psi ^{ - 2} - \psi ^{ - m - 1} \times \left( {1 + \psi ^{ - m}} \right)^{\frac{{1 - m}}{m}}{,}$$
(8)
$$S_m = \frac{{\partial Z}}{{\partial m}} = - \left( {1 + \psi ^{ - m}} \right)^{\frac{1}{m}} \times \left( {\frac{1}{{m^2}} \times \ln \left( {1 + \psi ^{ - m}} \right) + \frac{1}{m} \times \psi ^{ - m} \times \frac{{\ln \psi }}{{1 + \psi ^{ - m}}}} \right){.}$$
(9)

### Sensitivities of R/P to ϕ (Sϕ) and m (S*m)

The sensitivities of R/P to ϕ (Sϕ) and m (S*m) can be derived from Eq. (2) as:

$$S_\phi = \frac{{\partial Z}}{{\partial \phi }} = \left( {1 + \phi ^m} \right)^{\frac{1}{m}} \times \frac{{\phi ^{m - 1}}}{{1 + \phi ^m}} - 1{,}$$
(10)
$$S^ \ast _m = \frac{{\partial Z}}{{\partial m}} = \left( {1 + \varphi ^m} \right)^{\frac{1}{m}} \times \left( { - \frac{1}{{m^2}} \times \ln \left( {1 + \varphi ^m} \right) + \frac{1}{m} \times \varphi ^m \times \frac{{\ln \varphi }}{{1 + \varphi ^m}}} \right){.}$$
(11)

## References

1. 1.

Gudmundsson, L., Greve, P. & Seneviratne, S. I. Correspondence: Flawed assumptions compromise water yield assessment. Nat. Commun. 8, https://doi.org/10.1038/Ncomms14795 (2017).

2. 2.

Zhou, G. Y. et al. Global pattern for the effect of climate and land cover on water yield. Nat. Commun. 6, https://doi.org/10.1038/Ncomms6918 (2015).

3. 3.

Fuh, B. On the calculation of the evaporation from land surface. Sci. Atmos. Sin. 5, 23–31 (1981).

## Acknowledgments

This study was supported by Natural Science Foundation of China (Project numbers: 41430529, 31770493, and 31770759) and GDASCX-0701.

## Author information

P.Z., Q.L., G.Z., X.W., and M.Z. designed the study. P.Z. and Q.L. drafted the manuscript. P.Z., Q.L., G.Z., X.W., M.Z., Z.L., X.C., X.L., Y.X., G.S., D.S., S.Z., L.H., and Y. S. contributed to the revision of the manuscript.

Correspondence to Mingfang Zhang.

## 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

Zhou, P., Li, Q., Zhou, G. et al. Reply to 'Flawed assumptions compromise water yield assessment'. Nat Commun 9, 4788 (2018) doi:10.1038/s41467-018-07065-5

• ### Quantifying the biophysical effects of forests on local air temperature using a novel three-layered land surface energy balance model

• Yongxian Su
• , Liyang Liu
• , Jianping Wu
• , Xiuzhi Chen
• , Jiali Shang
• , Philippe Ciais
• , Guoyi Zhou
• , Raffaele Lafortezza
• , Yingping Wang
• , Wenping Yuan
• , Yilong Wang
• , Hongou Zhang
• , Guangqing Huang
•  & Ningsheng Huang

Environment International (2019)


text:
# How do I show that given $p(A \cap B) > 0$, when $p(A \cap B)$ decreases, and $p( B)$ remains the same, $p(A)$ must decrease?

Just asking for a formal proof.

Assume $$p(A \cap B) > 0$$. And assume that there can be four possibilities: $$(A \cap B)$$, $$(A^{c} \cap B)$$, $$(A \cap B^{c})$$, $$(A^{c} \cap B^{c})$$. Let’s say $$p(A \cap B)$$ decreases. But $$p(B)$$ stays the same. How to I prove that $$p(A)$$ must decrease? Or is this wrong?

I’m thinking of using $$p(A \cap B) = p(A) + p(B) - p(A \cup B)$$. But I don’t know how to proceed.

Thanks.

• Seems like a good start. Are $a$ and $b$ independent? If not, this isnt necessarily true. For texxing variables, logical symbols, etc. try enclosing the following symbols inside of dollar signs like so: \$P(a \cap b)\$, \$P(a \cup b)\$, \$\lnot a \land b\$, \$a \land \sim b\$ – David Diaz Mar 14 at 22:14
• Welcome to Mathematics Stack Exchange! A quick tour will enhance your experience. Here are helpful tips to write a good question and write a good answer. For equations, please use MathJax. – dantopa Mar 14 at 22:16

The claim is true if $$A,B$$ are independent (both before and after the "decrease"), because then $$P(A\cap B) = P(A)P(B)$$. So holding $$P(B)$$ constant, $$P(A)$$ and $$P(A \cap B)$$ must decrease together, in fact by the same proportion.
The claim can be false if $$A,B$$ are dependent. Here is a counterexample:
• Before: $$P(A) = P(B) = 0.3, P(A \cap B) = 0.09$$ (they were independent)
• After: $$P(B) = 0.3, P(A) = 0.7, P(A \cap B) = 0$$ (now they become complementary: $$A = B^c$$)


text:
The ACT course is an educator-created,
tutor-tested ACT prep course designed to help you boost your ACT score.

### Students enroll every day to study, learn, and reach their goals.

Join them and start reaching your own goals:

### WHAT MAKES US DIFFERENT (BETTER)

• #### 0% Screen Scribble

Everyone is a visual learner, so many of our practice questions come with a video explanation. These aren't screen-scribbles or tripod-teachers. Each video is professionally animated and designed to help you truly understand the answer.

• #### Comprehensive Content

Focused practice questions mean that if you complete at least 75% of the course, you're hitting all the content areas you'll see on the test and should be able to tackle any question you encounter.

• #### 100% Your Pace

Because the course is online and self-guided, you can study for the ACT whenever and wherever you'd like.

### ENROLL FOR FREE

#### 3 MONTH ENROLLMENT DURATION

• #### $390.00 3 months of enrollment ### COURSE CURRICULUM • 01 ##### Welcome Show Content • Welcome to The Olive Book! • 02 ##### Discussion Show Content • Using the "Discussions" Feature • 03 ##### Math Show Content • Math Pre-Test • Question 1 • Question 2 • Question 3 • Question 4 • Question 5 • Question 6 • Question 7 • Question 8 • Question 9 • Question 10 • Question 11 • Question 12 • Question 13 • Question 14 • Question 15 • Question 16 • Question 17 • Question 18 • Question 19 • Question 20 • Question 21 • Question 22 • Question 23 • Question 24 • Question 25 • Question 26 • Question 27 • Question 28 • Question 29 • Question 30 • Question 31 • Question 32 • Question 33 • Question 34 • Question 35 • Question 36 • Question 37 • Question 38 • Question 39 • Question 40 • Question 41 • Question 42 • Question 43 • Question 44 • Question 45 • Question 46 • Question 47 • Question 48 • Question 49 • Question 50 • Question 51 • Question 52 • Question 53 • Question 54 • Question 55 • Question 56 • Question 57 • Question 58 • Question 59 • Question 60 • Question 61 • Question 62 • Question 63 • Question 64 • Question 65 • Question 66 • Question 67 • Question 68 • Question 69 • Question 70 • Question 71 • Question 72 • Question 73 • Question 74 • Question 75 • Question 76 • Question 77 • Question 78 • Question 79 • Question 80 • Question 81 • Question 82 • Question 83 • Question 84 • Question 85 • Question 86 • Question 87 • Question 88 • Question 89 • Question 90 • Question 91 • Question 92 • Question 93 • Question 94 • Question 95 • Question 96 • Question 97 • Question 98 • Question 99 • Question 100 • Question 101 • Question 102 • Question 103 • Question 104 • Question 105 • Question 106 • Question 107 • Question 108 • Question 109 • Question 110 • Question 111 • Question 112 • Question 113 • Question 114 • Question 115 • Question 116 • Question 117 • Question 118 • Question 119 • Question 120 • Question 121 • Math Progress Test • Question 122 • Question 123 • Question 124 • Question 125 • Question 126 • Question 127 • Question 128 • Question 129 • Question 130 • Question 131 • Question 132 • Question 133 • Question 134 • Question 135 • Question 136 • Question 137 • Question 138 • Question 139 • Question 140 • Question 141 • Question 142 • Question 143 • Question 144 • Question 145 • Question 146 • Question 147 • Question 148 • Question 149 • Question 150 • Question 151 • Question 152 • Question 153 • Question 154 • Question 155 • Question 156 • Question 157 • Question 158 • Question 159 • Question 160 • Question 161 • Question 162 • Question 163 • Question 164 • Question 165 • Question 166 • Question 167 • Question 168 • Question 169 • Question 170 • Question 171 • Question 172 • Question 173 • Question 174 • Question 175 • Question 176 • Question 177 • Question 178 • Question 179 • Question 180 • Question 181 • Question 182 • Question 183 • Question 184 • Question 185 • Question 186 • Question 187 • Question 188 • Question 189 • Question 190 • Question 191 • Question 192 • Question 193 • Question 194 • Question 195 • Question 196 • Question 197 • Question 198 • Question 199 • Question 200 • Question 201 • Question 202 • Question 203 • Question 204 • Question 205 • Question 206 • Question 207 • Question 208 • Question 209 • Question 210 • Question 211 • Question 212 • Question 213 • Question 214 • Question 215 • Question 216 • Question 217 • Question 218 • Question 219 • Question 220 • Question 221 • Question 222 • Question 223 • Question 224 • Question 225 • Question 226 • Question 227 • Question 228 • Question 229 • Question 230 • Question 231 • Question 232 • Question 233 • Question 234 • Question 235 • Question 236 • Question 237 • Question 238 • Question 239 • Question 240 • Question 241 • Question 242 • Question 243 • Question 244 • Question 245 • Question 246 • Question 247 • Question 248 • Question 249 • Question 250 • Question 251 • Question 252 • Question 253 • Question 254 • Question 255 • Question 256 • Question 257 • Question 258 • Question 259 • Question 260 • Question 261 • Question 262 • Question 263 • Question 264 • Question 265 • Question 266 • Question 267 • You've Got This • Math Post-Test • Math Tips: Strategies for ACT Math • Math Tips: ACT Math Formulas to Make it Easier • Math Tips: What Calculator To Use on the ACT • 04 ##### Calculator Tutorial Show Content • How to Use the Calculator (Part 1) • How to Use the Calculator (Part 2) • 05 ##### Math Practice Test Show Content • Before You Begin • ACT Math Practice Test -- Questions • ACT Practice Test -- Bubble Sheet • ACT Math Practice Test -- Answers and Explanations • 06 ##### English -- Passage I: Grammar Matters Show Content • English Pre-Test • Question 1 • Question 2 • Question 3 • Question 4 • Question 5 • Question 6 • Question 7 • Question 8 • Question 9 • Question 10 • Question 11 • Question 12 • Question 13 • Question 14 • Question 15 • Question 16 • Question 17 • Question 18 • Question 19 • Question 20 • Question 21 • Question 22 • 07 ##### English -- Passage II: Test Prep Payoff Show Content • Question 1 • Question 2 • Question 3 • Question 4 • Question 5 • Question 6 • Question 7 • Question 8 • Question 9 • Question 10 • Question 11 • Question 12 • Question 13 • Question 14 • Question 15 • 08 ##### English -- Passage III: Dreaming of Hawaii Show Content • Question 1 • Question 2 • Question 3 • Question 4 • Question 5 • Question 6 • Question 7 • Question 8 • Question 9 • Question 10 • Question 11 • Question 12 • Question 13 • Question 14 • Question 15 • English Post-Test • English Tips: Grammar Things to Remember • 09 ##### Reading Show Content • Reading Pre-Test • Reading In Your Head 0: Planting the Seed • Passage I (Prose): Twelve Years a Slave • Reading In Your Head I: Prose Pro-tip • Passage II (Natural Science): Aquatic Insects • Reading In Your Head II: Natural Science Pro-tip • Passage III (Social Science): Ford Pinto • Reading In Your Head III: Social Science Pro-tip • Passage IV (Humanities): Art Theory • Reading In Your Head IV: Humanities Pro-tip • Passage V (Paired Passage): Federalists • Reading In Your Head V: Paired Passage Pro-tip • Reading In Your Head VI: Recap of All Passages • Reading Post-Test • Reading Tips: Strategies for the ACT Reading Section • 10 ##### Reading Practice Test Show Content • Before You Begin • ACT Reading Practice Test -- Questions • ACT Practice Test -- Bubble Sheet • ACT Reading Practice Test -- Answers and Explanations • 11 ##### Science Show Content • Science Pre-Test • Passage I: Noise and Scores • Passage II: Conflicting Viewpoints • Passage III: Solubility • Passage IV: Solar Geometry • Passage V: Spirulina • Passage VI: Genetics • Passage VII: pH and Dissolved Oxygen • Passage VIII: Albedo • Passage IX: Day and Night Gym Temperatures • Science Post-Test • Science Tips: How The ACT Science Section is Structured • Science Tips: Science Literacy Study Guide • 12 ##### Essay Show Content • Explanation of the ACT Essay • Practice Prompt 1 • Evaluating Your Essay • ACT Writing Test Scoring Rubric • Sample Essay 1 • Practice Prompt 2 • Sample Essay 2 • Official ACT Essay Prompt • Essay Tips: Strategies for the ACT Essay • ACT Essay FAQs • 13 ##### Full Practice Tests Show Content • Free Official ACT Practice Test • 14 ##### Tips and Strategies Show Content • About the ACT: Timing, Structure, and How to Study • 4 Month Study Schedule • 3 Month Study Schedule • 2 Month Study Schedule • Week of ACT Checklist • Math Tips: Strategies for the ACT Math Section • Math Tips: ACT Math Formulas to Make it Easier • Math Tips: What Calculator To Use on the ACT • English Tips: Grammar Things to Remember • Reading Tips: Strategies for the ACT Reading Section • Science Tips: How The ACT Science Section is Structured • Science Tips: Science Literacy Study Guide • Essay Tips: Strategies for the ACT Essay • 15 ##### One Last Question Show Content • What did you think of the Olive Book course? ### TUITION #### 3 MONTH ENROLLMENT DURATION • ####$390.00

3 months of access


text:
# Introduction

This package implements a Bayesian changepoint detection algorithm of the flavor found in Stephens (1994). This approach assumes that given an average of several 2D crosscuts from a 3D bullet land scan which has had the global structure removed via the robust LOESS procedure Cleveland (1979), the land engraved area (LEA) is fairly flat, and the groove engraved area (GEA) possesses an approximately linear structure with a nonzero slope. An example of this basic structure is given below.

The idea behind the changepoint approach is that within either the left GEA, right GEA, or the LEA, the global structure is consistent and can either be described by a line with zero slope, a line with positive slope for the right GEA, or a line with negative slope for the left GEA. Finding the points where the GEAs and LEA meet is treated as a problem of model selection. That is, the best fitting statistical model, in terms of the magnitude of the likelihood, should be the one which assumes that the points at which the global structure changes align with where the GEAs and LEA meet. These points of global structural change are what we will call changepoints. Thus, our model will be defined in a piecewise fashion. In practice there are also complex additional patterns which may exist for a number of reasons, but this large scale structural assumption remains generally reasonable. The complex smaller scale patterns can be thought of as the dependence in the data after accounting for the global structure. Because of the nature of the model which we consider, it becomes necessary for computational reasons to perform a couple of data preprocessing steps. Specifically, we will scale the residuals from the robust LOESS procedure, and we will impute missing values. In the next section, we describe the model that we will use to identify changepoints, after which we will describe the estimation procedure which we use. Following that, we will describe the two main data preprocessing steps necessary for estimation to be computationally feasible.

# Model

Before introducing the model, we introduce some notation. First, let $$\{Y(x_i): i = 1,2, ..., n\}$$ denote the set of random variables representing the residuals from the robust LOESS procedure at the values $$x_i$$. For simplicity, also assume that $$x_1 < x_2 < ... < x_n$$. Also, let $$c_l$$ be the value of the left changepoint and $$c_r$$ be the value of the right changepoint. Here, the left changepoint is where the left GEA meets the LEA, and the right changepoint is where the right GEA meets the LEA. Also, denote the median centered $$x$$ values as $$x'_i = x_i - \tilde{x}$$ where $$\tilde{x}$$ is the median $$x$$ value. As mentioned in the previous paragraph, the complex small scale patterns, such as the striae, will be modeled through a covariance structure on the data that will be allowed to differ between each GEA and between the GEAs and LEA. We will construct the covariance matrices from the exponential covariance function $$K(t, t';\sigma, \ell) = \sigma^2 e^{-\frac{|t - t'|}{\ell}} = cov(Y(t), Y(t'))$$. The differences in covariance matrices for the GEAs and LEA will be reflected in the parameters $$\sigma$$ and $$\ell$$. The data model that we consider is then,

\begin{align} (Y(x_1), Y(x_2), ..., Y(x_{k_1})) &\sim N(\beta_{01}\mathbb{1} + \beta_{11} x'_{1:k_1}, \Sigma_1(\sigma_1, \ell_1)) \\ (Y(x_{k_1 + 1}), Y(x_{k_1 + 2}), ..., Y(x_{k_2})) &\sim N(0, \Sigma_2(\sigma_2, \ell_2)) \\ (Y(x_{k_2 + 1}), Y(x_{k_2 + 2}), ..., Y(x_n)) &\sim N(\beta_{02}\mathbb{1} + \beta_{12} x'_{k_2 + 1:n}, \Sigma_3(\sigma_3, \ell_3)), \end{align}

where $$x_{k_1} < c_l \leq x_{k_1 + 1}$$ and $$x_{k_2} < c_r \leq x_{k_2 + 1}$$ Here, $$x_{1:k}$$ denotes the column vector $$(x_1, x_2, ..., x_k)^\top$$, and $$\mathbb{1}$$ denotes the vector of ones. Independence is assumed between each of these three distributions for simplicity.

Thus the parameters that need to be estimated include the four mean parameters in the GEAs, the six covariance parameters (two for each of the three areas), and the two changepoint parameters, $$c_l$$ and $$c_r$$.

The above model encapsulates the essence of the approach. However, there are a few difficulties. The first difficulty is that there are not always two GEAs in a particular land. There may be one GEA, or the land may only consist of the LEA. Thus, the above model is actually conditional on there being two GEAs in the data. We also define models for when there is one GEA on the left, one GEA on the right, or no GEAs. The models are defined in an essentially identical way. Conditional on there being only one GEA, the left GEA model is defined as,

\begin{align} (Y(x_1), Y(x_2), ..., Y(x_{k})) &\sim N(\beta_{0}\mathbb{1} + \beta_{1} x'_{1:k}, \Sigma_1(\sigma_1, \ell_1)) \\ (Y(x_{k + 1}), Y(x_{k + 2}), ..., Y(x_{n})) &\sim N(0, \Sigma_2(\sigma_2, \ell_2)), \end{align}

and the right GEA model is defined as,

\begin{align} (Y(x_{1}), Y(x_{2}), ..., Y(x_{k})) &\sim N(0, \Sigma_1(\sigma_1, \ell_1)) \\ (Y(x_{k + 1}), Y(x_{k + 2}), ..., Y(x_n)) &\sim N(\beta_{0}\mathbb{1} + \beta_{1} x'_{k + 1:n} \Sigma_2(\sigma_2, \ell_2)). \end{align}

Finally, conditional on there being no GEAs in the data, the model is simply

\begin{align} (Y(x_{1}), Y(x_{2}), ..., Y(x_{n})) &\sim N(0, \Sigma(\sigma, \ell)). \end{align}

Thus, estimating the changepoint locations also involves selecting the most appropriate model. In order to avoid confusion, we have slightly abused notation and, for example, $$\Sigma_1(\sigma_1, \ell_1)$$ as it is estimated in the two changepoint model is not the same as $$\Sigma_1(\sigma_1, \ell_1)$$ from either of the one changepoint models, and $$\Sigma_1(\sigma_1, \ell_1)$$ is also not the same between the two one changepoint models. As another example, $$\beta_0$$ is not the same between each of the one changepoint models. So, to be clear, duplication of notation in different models is not meant to imply that those parameters are shared between models.

Ultimately, these above four models are each individually fitted, and each model above is given a prior. From there, we do model selection in the formal Bayesian way, selecting the most probable model. Simultaneously with selecting the most probable model, we also use the maximum a posteriori estimator for the changepoint locations.

In order to complete a Bayesian model specification, we need priors on each of the parameters in each model as well as each model. We will assume independence between each parameter a priori. For each length scale $$\ell$$, we will assume $$\ell \sim \text{Gamma}(3,5)$$. For each standard deviation, we will assume $$\sigma \sim \text{Half-Normal}^{+}(0,1)$$, where $$\text{Half-Normal}^{+}(\cdot,\cdot)$$ is notation for the normal distribution restricted to the positive real numbers. For intercept parameters, $$\beta_{01}, \beta_{02}, \beta_0 \sim N(0, 10)$$. For the slope parameters, the preceding trend deviates slightly. For any slope that corresponds to the left GEA, $$\beta_1$$ or $$\beta_{01}$$, we will assume that the slope can not be positive. That is, $$\beta_1, \beta_{01} \sim \text{Half-Normal}^{-}(0,10)$$, where $$\text{Half-Normal}^{-}(\cdot, \cdot)$$ is notation for the normal distribution restricted to the negative real numbers. Contrastingly, for any slope that corresponds to the right GEA, $$\beta_1$$ or $$\beta_{02}$$, we will assume that the slope can not be negative. That is, $$\beta_1, \beta_{01} \sim \text{Half-Normal}^{+}(0,10)$$. For the changepoint locations, we assume a uniform prior $$\pi(c_l, c_r) \propto I(a < c_l < c_r - \gamma < b - \gamma)$$. Here, $$a$$ and $$b$$ are some values close to the edges of the data. How close those values are to the edges is a parameter that is set manually. Further, we include another hyperparameter, $$\gamma$$, which can be set so that the changepoints are not allowed to be too close to each other. This is also a parameter that is set manually. Lastly, we assume a uniform prior over all four models.

# Estimation

As was noted in Stephens (1994), for any model including a changepoint, the likelihood is not a smooth function of the changepoint location. This is because, holding all other parameters fixed, shifting the changepoint value will result in zero change to the likelihood until it crosses the nearest point to the right or left, at which point the likelihood makes a jump. This makes maximum likelihood estimation in the standard way infeasible, but Bayesian estimation can be done in a fairly straightforward way via Gibbs sampling. The basic idea is that, for each model, we can construct a two step Gibbs sampler. In step 1 we sample from the posterior distribution of the mean and covariance parameters given the changepoint locations, and in step 2 we sample from the changepoint locations given the mean and covariance parameters. Because of the non-conjugacy in our model, we perform both sampling steps using a random walk Metropolis-Hastings (RWMH) step with Gaussian proposals. For details on Gibbs sampling and the Metropolis-Hastings algorithm see Gelman et al. (2013). It is also worth mentioning that the zero changepoint model does not require Gibbs sampling at all, and we perform estimation there using a RWMH algorithm.

As a practical note, it turns out that the posterior distribution is almost always multimodal, and it can happen that the sampler gets stuck in a suboptimal mode for a large number of iterations. It is also the case that the suboptimal modes need not even be close to the groove locations. It has, however, been our experience that the optimal mode corresponds well to the actual groove locations, which are often somewhat close to the edges of the data. With this in mind, starting values and the RWMH proposal variances play a very important role in the success of the sampling algorithm. Fortunately, it seems to be the case that by setting the initial changepoint values close to the edges of the data and making the proposal variance small (around 100 seems to work well) allows the sampler to wander inwards, and even with a modest number of iterations (say 5000), typically pass through the largest mode corresponding to the groove locations. This is not always the case, and it is possible that increasing the number of iterations produces better results.

The sampling functions were originally written with the intention of tuning the proposal variances to potentially accelerate convergence, and thus several warmup iterations are required for this purpose. This turns out to be a bad idea in this context for two reasons. The first reason is that the warmup iterations allow the sampler to wander past the global modes and get stuck in suboptimal modes far from the groove locations, from which the sampler may or may not find its way back to the optimal modes in just a few thousand iterations. Secondly, if the sampler does wander past the optimal modes, which are usually on the edges of the data, the tuned proposal variance can be quite large. The large proposal variance might not be a huge problem if it weren’t for the fact that the width of the modes are almost always quite small. This means that it can take a very, very long time for the sampler to move from a suboptimal mode to the global mode. In order to mitigate this problem, we are currently setting the number of warmup iterations to be relatively small (somewhere in 100 to 500 seems to work well). In future iterations of these functions, we will remove the necessity of having warmup iterations or tuned proposal variances.

We now provide the two basic steps of the Gibbs sampler for the two changepoint case. The algorithms to sample from the other three models are omitted, and are nearly identical except for the a smaller number of parameters need to be sampled. Denote collection of mean and covariance parameters for the left GEA as $$\theta_1$$, the LEA as $$\theta_2$$, and the right GEA as $$\theta_3$$. Then, at iteration $$t$$ after warmup

1. given changepoint locations $$(c_l^{(t - 1)}, c_r^{(t - 1)})$$, sample $$(\theta_1^{(t)}, \theta_2^{(t)}, \theta_3^{(t)})$$ using independent RWMH steps for each $$\theta_i$$
2. given $$(\theta_1^{(t)}, \theta_2^{(t)}, \theta_3^{(t)})$$, sample $$(c_l^{(t)}, c_r^{(t)})$$ using a single RWMH step.

Initially, the Metropolis proposal variance for each $$\theta_i$$ is diagonal with diagonal elements all equal to $$1/2$$. The proposal variance for $$(c_l, c_r)$$ is initially set to be diagonal with elements equal to $$10^2$$. Note that because of the currently necessary warmup iterations, the variances after warmup for each $$\theta_i$$ becomes $$\frac{2.4^2}{d}\hat{Var}(\theta_i^{(1:w)}) + \text{diag}(0.1)$$, where $$d$$ is the dimension of $$\theta_i$$ (which is not constant between GEAs and LEA), and $$\hat{Var}(\theta_i^{(1:w)})$$ is the estimated variance covariance matrix from the $$w$$ warmup iterations. Note that the addition of a diagonal matrix with entries $$0.1$$ is to avoid the case when most or all warmup iterations have the same value. Similarly, the proposal variance for $$(c_l, c_r)$$ after warmup becomes $$\frac{2.4^2}{2}\hat{Var}((c_l,c_r)^{(1:w)}) + \text{diag}(1)$$.

After running the MCMC for each model, parameter estimates and the most likely model are jointly chosen according to the largest joint posterior value. That is, we arrive at estimates $$(\hat{\theta}, \hat{M}) = \underset{(\theta, M)}{\operatorname{argmax}}{\log(p(\theta, M | Y))}$$, where $$M$$ is the random variable associated with the choice of model, $$\theta$$ is the associated parameter vector for the appropriate model, and $$Y$$ is all of the available data.

# Data Preprocessing

Before running the analysis described above, we first perform two data preprocessing steps. The first step is to scale the residuals from the robust loess procedure by the standard deviation calculated from the entire set of residuals. The reason for this is simply to make priors for standard deviation and slope parameters easier to specify. For example, ensuring that the residuals are scaled to have standard deviation one means that the standard deviation parameters in our model should also be close to one. This scaling also ensures that slopes values are not very large.

The second preprocessing step is a bit more involved. In order to enable the algorithm to run reasonably fast, we need to take advantage of the sparse precision matrix structure that is induced by the covariance function we chose. Indeed, this was the reason for choosing this covariance function in the first place. Unfortunately, it is challenging to do this unless the observations are evenly spaced in the domain. In our case, this would be true if there were no missing values. In order to remedy this problem, we impute the missing data, but only in the case that there exist nonmissing observations outside of the missing values. In the case that the missing values exist on the edges of the data, we simply do not consider those domain values in the model.

We perform the imputation by treating the observations as coming from an unknown function, and infer the missing values from the known function values. In order to do this, we model the data with a Gaussian process and the squared exponential covariance function. That is, we suppose that

$Y(t) \sim \mathcal{GP}(0, K(t,t';\sigma^2, \ell)),$

where now $$K(t,t';\sigma^2, \ell) = \sigma^2 e^{-(t - t')^2/(2\ell^2)}$$ is the squared exponential covariance function. We emphasize for clarity that this is a different covariance function than we use in the changepoint model. The main reason for this is that in imputing values, it seems desirable to allow dependencies beyond immediately neighboring points to influence predictions as the function that we are trying to predict generally has a smooth global structure.

In the changepoint algorithm, we do allow the user to choose to estimate the Gaussian process parameters via maximum likelihood, which just uses the optim() function. However, we do not use the full set of data to estimate the parameters, as there is still too much to do estimation reasonably quickly. For this reason, we take a subset of the data chosen by selecting every 20-th point starting at the first observed location. This still provides plenty of data to estimate the covariance parameters well. However, we also supply default covariance parameter values so that estimation of these parameters is not necessary. These values were estimated from a representative bullet land, and we suspect that they will almost always work sufficiently well in practice.

When we impute the missing values, we use the estimated or default covariance parameter values and all of the observed data and compute the conditional mean of the missing values. To be clear, denote the distribution of the observed and missing data as

$(Y,Y^*)^\top \sim N\left( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} \Sigma_{yy} & \Sigma_{yy^*} \\ \Sigma_{y^*y} & \Sigma_{y^*y^*}\end{bmatrix} \right).$

Here, $$Y$$ is observed data and $$Y^*$$ is the missing data, and the covariance matrix above is constructed from the squared exponential covariance function and the estimated or default covariance parameter values. We then use normal distribution theory to calculate the imputed values

$E(Y^*|Y = y) = \Sigma_{y^*y} \Sigma_{yy}^{-1}y$.

# Example

In this section, we will go through an example by individually using each necessary function as well as by using the wrapper function designed to perform all of the necessary steps with one function call. Data from this example comes from the the James Hamby study (Hamby, Brundage, and Thorpe 2009). We downsample the data to make the example run faster, but this would not be preferred in practice.

First, we read in the data. These data are an example of an average of several 2D cross cuts from a 3D bullet land scan in the Hamby 44 data set. We can tell both from the plot and the head(raw_data) call that there are missing values in the data.

library(bulletcp)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#>     filter, lag
#> The following objects are masked from 'package:base':
#>
#>     intersect, setdiff, setequal, union

data("example_data")
#>       x value
#> 1 0.000    NA
#> 2 0.645    NA
#> 3 1.290    NA
#> 4 1.935    NA
#> 5 2.580    NA
#> 6 3.225    NA
raw_data <- raw_data[seq(from = 1, to = nrow(raw_data), by = 30),]

ggplot(data = raw_data) +
geom_point(aes(x = x, y = value)) +
theme_bw() +
ylab("Height") +
xlab("Width")
#> Warning: Removed 11 rows containing missing values (geom_point).

## Removing global structure

The data shown above still contains the global bullet curvature. The changepoint algorithm requires removal of this global structure before imputing the missing data. Functionality to do this is provided in the bulletxtrctr package (Hofmann, Vanderplas, and Krishnan 2018). Code to remove this structure is shown below.

# set the minimum y-value to zero
check_min <- min(raw_data$value[!is.na(raw_data$value)])
raw_data <- dplyr::mutate(raw_data, value_std = value - check_min)

# remove global structure
rlo_fit <- robust_loess_fit(cc = raw_data, iter = 20)
raw_data$rlo_pred <- predict(rlo_fit, newdata = raw_data) raw_data$rlo_resid <- raw_data$value_std - raw_data$rlo_pred

# define new data frame without the global structure
data <- data.frame("x" = raw_data$x, "y" = raw_data$rlo_resid)

# plot the new data
ggplot(data = data) +
geom_point(aes(x = x, y = y)) +
theme_bw() +
ylab("Height") +
xlab("Width")
#> Warning: Removed 11 rows containing missing values (geom_point).

## Imputing missing values

Now that the global structure has been removed from the land, we can proceed to impute missing values. We will work with scaled y-values, and so the first step is to perform this scaling.

# make range of x data equal to the range of non NA points and scale the y-values
d <- data.frame("x" = data$x, "y" = scale(data$y))
max_x_notNA <- max(d$x[!is.na(d$y)])
min_x_notNA <- min(d$x[!is.na(d$y)])
d <- d[d$x >= min_x_notNA & d$x <= max_x_notNA,]

Once the data have been scaled, we take a subset of them by selecting every 20-th point to use when estimating the Gaussian process parameters. Then, we estimate these parameters using the mlgp() function in this package. The optimization routine is done on the log scale so as to remove boundary constraints for the Gaussian process parameters, which are both constrained to be positive. The first parameter value returned is the log of the process standard deviation, and the second returned parameters value is the log of the length scale. The default values used if one does not choose to estimate the Gaussian process parameters are $$0.8$$ and $$15$$ for the standard deviation and length scale, respectively.

# subset data
temp_d <- d[seq(from = 1, to = nrow(d),by = 1),]

# remove missing data
temp_d <- temp_d[complete.cases(temp_d),]

# estimate the MLE's
mles <- mlgp(y = temp_d$y, x = temp_d$x)

# name MLE's impute_par
impute_par <- exp(mles$par) impute_par #> [1] 0.9671305 22.2714492 Now that we have parameter estimates, we can use the function myimpute() to impute missing values. full_data <- imputeGP(y = d$y, x = d$x, sigma = impute_par[1], l = impute_par[2]) head(full_data) #> x y #> 1 38.70 7.107245403 #> 2 58.05 5.181357035 #> 3 77.40 1.568790923 #> 4 96.75 1.831877298 #> 5 116.10 -0.000138855 #> 6 135.45 -0.337581864 ggplot(data = full_data) + geom_point(aes(x = x, y = y)) + theme_bw() ## Estimating changepoint locations Now that we have imputed missing values, we can proceed to estimate the number and location of any grooves. In order to do this, we can can either call the function variable_cp_gibbs_v2(), which, in turn, calls functions for running separate Gibbs sampling algorithms in the 0, 1, or 2 changepoint cases, or we can call each of these functions individually. To clearly demonstrate what is going on, we will call each of the necessary functions directly. The only things that we really care about are the estimated groove locations, as well as the maximum value of the log posterior density (which occurs at the estimated groove locations). The latter will enable us to choose the number and location of any possible grooves while the former actually provides the estimated groove locations in each model. ### Model with no GEAs First, we will run a Gibbs sampler estimating model parameters in the case that there are zero changepoints. In this case, there should be no GEAs in the data. This is clearly not true in our example, but in some cases it will be. The function which runs the MCMC sampler in this case is named cp0_gibbs(). It may be worth noting that, although the function name would seem to imply that the algorithm is a Gibbs sampler, in the case of zero changepoints, it is actually just a RWMH algorithm. The function was named for the sake of consistency across sampling functions, which are always Gibbs samplers except in this special case. Argument descriptions for this function can be found in the relevant function documentation. In the beginning of the following code chunk, a function named lognormal_ou_pdf() is defined. This function defined the log of the pdf of a normal distribution with the covariance structure in the models described above. This function is necessary as it avoids matrix operations, taking advantage of the particular covariance structure we have defined. It is not critical to understand this function beyond the fact that it is needed in every MCMC sampler we employ. # log pdf of a normal distribution with the type of covariance matrix described in the model lognormal_ou_pdf <- function(x, mu, sigma, l) { n <- length(x) rho <- exp(-1/l) return(-n/2 * log(2 * pi) - n * log(sigma) - ((n - 1)/2) * log(1 - rho^2) - 1/2 * 1/(sigma^2 * (1 - rho^2)) * ((x[1] - mu[1])^2 + (x[n] - mu[n])^2 + (1 + rho^2) * sum((x[2:(n-1)] - mu[2:(n-1)])^2) - 2 * rho * sum((x[1:(n-1)] - mu[1:(n-1)]) * (x[2:n] - mu[2:n])))) } # define starting values start.vals <- list("sigma" = c(1), "l" = c(10)) # proposal variance for the MH step prop_var <- diag(c(1/2,1/2)) # run the RWMH sampler for the zero changepoint model set.seed(1111) m0cp <- runmcmc_cp0(data = full_data, iter = 1000, start.vals = start.vals, prop_var = prop_var, warmup = 500, verbose = FALSE) # print maximum value of the (conditional) log posterior up to the normalizing constant max(m0cp$lpost)
#> [1] NA

### Model with one GEA

Now we proceed to the slightly more involved cases where the model assumes that there are one or more changepoints in the data. If we assume that there is one changepoint, we must model the case when the GEA is on the left side of the land as well as on the right side of the land. The functions that we will use to do this are named runmcmc_cp1_left() and runmcmc_cp1_right(), corresponding to the models with the location of the the GEA assumed to be on the left and right side of the land, respectively. We must specify more arguments to these functions than we did in the previous case. These arguments mostly correspond to parameters for the proposal variances for sampling the parameters on both the GEA and the LEA, as well as starting values for each of these sets of parameters. However, we must also specify the tol_edge argument which forces the potential changepoint to be a minimum distance from the edge of the data.

# define starting values for the changepoints
cp_start_left <- min(full_data$x) + 60 cp_start_right <- max(full_data$x) - 60

# define list of starting values for both the left and right changepoint models
start.vals <- list("left" = list("sigma" = c(1,1),
"l" = c(10,10),
"cp" = c(cp_start_left),
"beta" = c(-1),
"intercept" = c(0)),
"right" = list("sigma" = c(1,1),
"l" = c(10,10),
"cp" = c(cp_start_right),
"beta" = c(1),
"intercept" = c(0)))

# list of starting values for each of the two MH steps (not sampling the changepoint) for both the left and right changepoint models
prop_var <- list("left" = list(diag(c(1/2,1/2,1/2,1/2)),
diag(c(1/2,1/2))),
"right" = list(diag(c(1/2,1/2)),
diag(c(1/2,1/2,1/2, 1/2))))

# define the proposal variance for the RWMH step sampling the changepoint
cp_prop_var <- 10^2

# run Gibbs MCMC for the left GEA model
set.seed(1111)
m1cp_left <- runmcmc_cp1_left(data = full_data, iter = 1000, warmup = 500,
start.vals = start.vals$left, prop_var = prop_var$left,
cp_prop_var = cp_prop_var,
verbose = FALSE, tol_edge = 50)

# run Gibbs MCMC for the right GEA model
set.seed(1111)
m1cp_right <- runmcmc_cp1_right(data = full_data, iter = 1000, warmup = 500,
start.vals = start.vals$right, prop_var = prop_var$right,
cp_prop_var = cp_prop_var,
verbose = FALSE, tol_edge = 50)

# Find MAP estimates under both models of the changepoint location
map_cp1_left <- as.numeric(m1cp_left$parameters$cp)[which.max(m1cp_left$lpost)] map_cp1_right <- as.numeric(m1cp_right$parameters$cp)[which.max(m1cp_right$lpost)]

# print maximum value of (conditional) log posterior up to a normalizing constant
max(m1cp_left$lpost) #> [1] -39.2708 max(m1cp_right$lpost)
#> [1] -38.94021

# Plot the data and the changepoint MAPs
ggplot(data = full_data) +
geom_point(aes(x = x, y = y)) +
geom_vline(aes(xintercept = map_cp1_left)) +
geom_vline(aes(xintercept = map_cp1_right)) +
theme_bw()

In the plot above, the vertical black lines correspond to the MAP estimates of the changepoint locations for each model. Also, it is clear that the left changepoint model is a better fit to the data given the fact that the maximum log posterior value (up to a constant) is larger than it is for the right changepoint model. However, both of these models appear more likely than the zero changepoint model, which has a much lower maximum log posterior value.

Note that the function runmcmc_cp1() is just a wrapper for the runmcmc_cp1_left() and runmcmc_cp1_right() functions above. To use this function, starting values and proposal variances for each model are input separately. As shown below, the MAP estimates for the changepoints are very close to the estimates found from running each analysis individually as above. They will never be identical without setting a seed, as the MCMC algorithm is inherently random.

# run both Gibbs MCMC algorithms for each of the single changepoint model
set.seed(1111)
m1cp <- runmcmc_cp1(data = full_data, iter = 1000,
start.vals.left = start.vals$left, start.vals.right = start.vals$right,
prop_var_left = prop_var$left, prop_var_right = prop_var$right,
cp_prop_var = cp_prop_var,
tol_edge = 50,
warmup = 500, verbose = FALSE)

# Find MAP estimates under both models of the changepoint location
map_cp1_left2 <- as.numeric(m1cp$left_parameters$cp)[which.max(m1cp$lpost$left)]
map_cp1_right2 <- as.numeric(m1cp$right_parameters$cp)[which.max(m1cp$lpost$right)]

map_cp1_left2
#> [1] 153.8055
map_cp1_right2
#> [1] 165.0584

### Model with two GEAs

Because assuming a model structure with a single GEA involves fitting two models, fitting a model that assumes there are two GEAs is actually a bit simpler. Though it requires proposal variances and starting values for both GEAs as well as the LEA, it requires only performing one MCMC. The function that performs this sampling is named changepoint_gibbs_2cp_v2(), and its use is shown in the code block below. The one additional argument specific to this function is the named tol_cp. This is similar to the tol_edge argument except it specifies how close the changepoints can be to each other, not to the edge of the data. The default value for this argument is 1000, as this range is well below the range of the data.

# define starting values
start.vals <- list("sigma" = c(1,1,1),
"l" = c(10,10,10),
"cp" = c(cp_start_left, cp_start_right),
"beta" = c(-2,2),
"intercept" = c(0,0))

# define proposal variances (not for changepoints)
prop_var <- list(diag(c(1/2,1/2,1/2,1/2)),
diag(c(1/2,1/2)),
diag(c(1/2,1/2,1/2,1/2)))

# define proposal variance for changepoints
cp_prop_var <- diag(c(10^2, 10^2))

# run Gibbs MCMC
set.seed(1111)
m2cp <- runmcmc_cp2(data = full_data,
iter = 1000,
start.vals = start.vals,
prop_var = prop_var,
cp_prop_var = cp_prop_var,
tol_edge = 50,
tol_cp = 1000,
warmup = 500,
verbose = FALSE)

# Find MAP estimates of changepoint locations
map_cp2 <- m2cp$parameters$cp[which.max(m2cp$lpost),] # print maximum value of (conditional) log posterior up to a normalizing constant max(m2cp$lpost)
#> [1] -35.88717

# Plot the data and the changepoint MAPs
ggplot(data = full_data) +
geom_point(aes(x = x, y = y)) +
geom_vline(aes(xintercept = map_cp2[1])) +
geom_vline(aes(xintercept = map_cp2[2])) +
theme_bw()

This model has the highest maximum log posterior value, and thus will be chosen as the appropriate model for the data. Furthermore, the estimated changepoints appear to be very close to the true locations where the GEAs and the LEA meet.

### Running all MCMCs with one function

Of course, everything done above can be done much more simply by using one of the three wrapper functions designed to make this entire process simple. There are wrapper functions which differ only in terms of the preprocessing steps that they perform and the number of objects which they return. All of the functions run all of the MCMCs and estimate groove locations. We will now show how each of these functions is designed to be used in practice. The only additional argument that is required by these functions is the prior on the number of grooves named prior_numcp. This argument is a four element vector that must sum to 1. The elements of the prior vector are ordered as such (LEA only model, left GEA only model, right GEA only model, two GEA model). In practice it seems that the differences in the log posteriors (conditional on the model) are large enough that the model prior would need to be extraordinarily strong to affect the final results.

First, we will start with the function which does the least data preprocessing. This function is named runmcmc_cpall(). It requires that all data preprocessing be done beforehand, and it only runs the MCMCs and estimates the groove locations. We also need to specify lists of starting values and proposal variances. Each list element corresponds to a given model, so specification of starting values and proposal variances is a little cumbersome, but default starting values seem to work reasonably well on cases that have been tested so far. The list elements need to be named as they are in the example code below.

# changepoint starting values
cp_sval_left <- min(full_data$x) + 60 cp_sval_right <- max(full_data$x) - 60

# list of starting values for each model
start.vals <- list("cp2" = list("sigma" = c(1,1,1),
"l" = c(10,10,10),
"cp" = c(cp_sval_left,cp_sval_right),
"beta" = c(-2,2),
"intercept" = c(0,0)),
"cp1" = list("left" = list(
"sigma" = c(1,1),
"l" = c(10,10),
"cp" = c(cp_sval_left),
"beta" = c(-1),
"intercept" = c(0)),
"right" = list(
"sigma" = c(1,1),
"l" = c(10,10),
"cp" = c(cp_sval_right),
"beta" = c(1),
"intercept" = c(0))),
"cp0" = list("sigma" = c(1), "l" = c(10)))

# proposal variances for each model
prop_var <- list("cp2" = list(diag(c(1/2,1/2,1/2,1/2)),
diag(c(1/2,1/2)),
diag(c(1/2,1/2,1/2,1/2))),
"cp1" = list("left" = list(diag(c(1/2,1/2,1/2,1/2)),
diag(c(1/2,1/2))),
"right" = list(diag(c(1/2,1/2)),
diag(c(1/2,1/2,1/2, 1/2)))),
"cp0" = diag(c(1/2,1/2)))

# changepoint proposal variances
cp_prop_var <- list("cp2" = diag(c(10^2, 10^2)),
"cp1" = 10^2)

# prior on the number of changepoints
prior_numcp <- rep(1/4, times = 4)

set.seed(1111)
cp_gibbs <- runmcmc_cpall(data = full_data,
start.vals = start.vals,
prop_var = prop_var,
cp_prop_var = cp_prop_var,
verbose = FALSE,
tol_edge = 50,
tol_cp = 1000,
iter = 1000,
warmup = 500,
prior_numcp = prior_numcp)

# print the estimated maximum log posteriors
cp_gibbs$max_lpost #>$cp0
#> [1] NA
#>
#> $cp1_left #> [1] -40.48337 #> #>$cp1_right
#> [1] -40.69286
#>
#> $cp2 #> [1] -37.27347 # Estimated groove locations in each model cp_gibbs$cp_map
#> $2cp #> [1] 110.4993 1993.0634 #> #>$1cp
#> $1cp$left
#> [1] 148.1259
#>
#> $1cp$right
#> [1] 127.2922

Now we introduce the usage of a function which, once global structure has been removed, will scale the data and impute missing values. The function is named detect_cp(). Notice that the object passed to the data argument is the unscaled data that still has missing values but has the global structure removed. Here, we have also set est_impute_par = TRUE indicating that we would like to estimate the imputation parameters. If this argument is set to FALSE, then the argument impute_par will need to be set with user supplied covariance parameters.

# Impute missing data, run MCMCs, and estimate MAPs
set.seed(1111)
cp_gibbs2 <- detect_cp(data = data,
start.vals = start.vals,
prop_var = prop_var,
cp_prop_var = cp_prop_var,
verbose = FALSE,
tol_edge = 50,
tol_cp = 1000,
iter = 1000,
warmup = 500,
prior_numcp = prior_numcp,
est_impute_par = TRUE)

# print the estimated log posteriors
cp_gibbs2$changepoint_results$max_lpost
#> $cp0 #> [1] NA #> #>$cp1_left
#> [1] -34.20846
#>
#> $cp1_right #> [1] -35.34848 #> #>$cp2
#> [1] -41.51201

# Estimated groove locations in each model
cp_gibbs2$changepoint_results$cp_map
#> $2cp #> [1] 106.5513 1991.7210 #> #>$1cp
#> $1cp$left
#> [1] 126.6153
#>
#> $1cp$right
#> [1] 184.6226

# The groove locations automatically chosen based on the MAP criteria
cp_gibbs2$grooves #> [1] 126.6153 2051.1000 Finally, we can call a function that also removes global structure, and thus can take in the raw crosscut data. This function is called get_grooves_bcp. The necessary arguments of this function follow a slightly different format than with the other functions so that it conforms to similar functions existing in the bulletxtrctr package. This function ultimately just removes the global structure and then calls detect_cp(). As such, the user can supply arguments to be passed to detect_cp, but it is intended that this not usually be necessary. The main reason why one may want to do this is to increase the number of iterations that each MCMC is run. The longer the MCMC is run, the more likely that the models can escape any local modes. For this reason, we demonstrate how a specific user supplied argument for the number of iterations can be fed to get_grooves_bcp. # Estimate the groove locations by supplying additional arguments cp_gibbs3 <- get_grooves_bcp(x = raw_data$x, value = raw_data$value, adjust = 10, iter = 1000) # Estimated groove locations cp_gibbs3$groove
#> [1]  108.8436 2041.1000

# Acknowledgements

We would like to thank the efforts of CSAFE and the Roy J Carver High Resolution Microscopy lab in scanning the Hamby 44 set and providing the scans to us.

# References

Cleveland, William S. 1979. “Robust Locally Weighted Regression and Smoothing Scatterplots.” Journal of the American Statistical Association 74 (368): 829–36.

Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. 2013. Bayesian Data Analysis. Chapman; Hall.

Hamby, James E., David J. Brundage, and James W. Thorpe. 2009. “The Identification of Bullets Fired from 10 Consecutively Rifled 9mm Ruger Pistol Barrels: A Research Project Involving 507 Participants from 20 Countries.” AFTE Journal 41 (2): 99–110.

Hofmann, Heike, Susan Vanderplas, and Ganesh Krishnan. 2018. Bulletxtrctr: Automatic Matching of Bullet Striae. https://heike.github.io/bulletxtrctr/.

Stephens, David A. 1994. “Bayesian Retrospective Multiple-Changepoint Identification.” Journal of the Royal Statistical Society. Series C (Applied Statistics) 43 (1): 159–78.


text:
## Stream: new members

### Topic: Application type mismatch

#### Vincent Beffara (Nov 22 2019 at 13:12):

Hi,

Still trying to figure out how lean works ... so I have a function defined under the condition that some predicate holds on the argument, something like

def good (x : whatever) : Prop := something
def f (x : whatever) (h : good x) := blah

and I would like to rewrite the argument x in f x h, but I get the error [check] application type mismatch at ... stating, if I understand the message correctly, that after rewriting the x to some y, then we would get f y h with h : good x instead of good y. Certainly doing the same rewrite inside of h would make it fit, but if I understand it right, rw only acts at one location at a time.

What is the right way to do this rewrite? Is that what simp only is for? Is that kind of issue the reason why 1/0=0?

#### Vincent Beffara (Nov 22 2019 at 13:13):

OK, more explicit example:

def good (n : ℕ) : Prop := 3 ∣ n
def f (n : ℕ) (h : good n) := n
lemma toto {n : ℕ} : n = n + 1 - 1 := by { norm_num }
example {n : ℕ} (h : good n) : 3 ∣ f n h :=
begin
rw (@toto n),
sorry
end

#### Johan Commelin (Nov 22 2019 at 13:24):

Is that kind of issue the reason why 1/0=0?

Yes, this is certainly one of the reasons

#### Rob Lewis (Nov 22 2019 at 13:24):

This example will work if you revert h before the rewrite. Then both occurrences of n will get rewritten at the same time. But in general, this will be frustrating, and yes it is the reason why 1/0=0.

#### Rob Lewis (Nov 22 2019 at 13:24):

In mathlib we generally prefer to "totalize" f, that is, to give it a default value when good n doesn't hold. Something like

#### Johan Commelin (Nov 22 2019 at 13:24):

@Vincent Beffara It turns out that it is often easier to define your function in such a way that it garbage whenever the input is garbage, instead of defining a partial function

#### Rob Lewis (Nov 22 2019 at 13:24):

@[derive decidable] def good (n : ℕ) : Prop := 3 ∣ n

def f (n : ℕ) : ℕ := if h : good n then n else default _

lemma f_val {n : ℕ} (h : good n) : f n = n := by simp [f, h]

lemma toto {n : ℕ} : n = n + 1 - 1 := by { norm_num }
example {n : ℕ} (h : good n) : 3 ∣ f n :=
begin
rw [(@toto n)] at h ⊢,
sorry
end

#### Johan Commelin (Nov 22 2019 at 13:26):

You only need to check good x on theorems that need it, but not on the function itself.
If you can make the "output garbage" somewhat useful even if the input is garbage, that's even better, because it will allow you to drop some checks of good x. That's the reason why 1/0 = 0 and not 1/0 = 37 even though Kevin would have liked that.

#### Kevin Buzzard (Nov 22 2019 at 13:28):

Here's another trick which is sometimes helpful:

import tactic.norm_num

def good (n : ℕ) : Prop := 3 ∣ n
def f (n : ℕ) (h : good n) := n
lemma toto {n : ℕ} : n = n + 1 - 1 := by { norm_num }
example {n : ℕ} (h : good n) : 3 ∣ f n h :=
begin
have : 3 ∣ f (n + 1 - 1) h, sorry,
convert this,
end

#### Kevin Buzzard (Nov 22 2019 at 13:29):

toto is true by rfl but had this not been the case, the goal would have changed to n = n + 1 - 1. Possibly.

#### Vincent Beffara (Nov 22 2019 at 16:16):

This example will work if you revert h before the rewrite. Then both occurrences of n will get rewritten at the same time. But in general, this will be frustrating, and yes it is the reason why 1/0=0.

Ah, that is nice to know. So how can I tell whether rw will rewrite the first occurrence of n or all of them at once?

And a related question, if the h is not part of my environment but a field in some structure somewhere (and the n is in another field of the same structure), does all that mean that I will have to deconstruct the whole structure before I can apply rw then?

#### Kevin Buzzard (Nov 22 2019 at 16:17):

You only need to check good x on theorems that need it, but not on the function itself.
If you can make the "output garbage" somewhat useful even if the input is garbage, that's even better, because it will allow you to drop some checks of good x. That's the reason why 1/0 = 0 and not 1/0 = 37 even though Kevin would have liked that.

It's morally dishonest to prove theorems about junk.

#### Kevin Buzzard (Nov 22 2019 at 16:18):

You can target which ns rw will rewrite, it's somewhat arcane but it's possible.

#### Vincent Beffara (Nov 22 2019 at 16:19):

Is it arcane even if I want to target "all the ns even those hidden in assumptions that show up as _"?

#### Kevin Buzzard (Nov 22 2019 at 16:20):

I just usually proceed by trial and error :-/ I try rw, erw, simp only and then I ask here :D

#### Vincent Beffara (Nov 22 2019 at 16:21):

Feels a bit brittle to navigate through the expression and tag all the occurrences one by one ...

#### Patrick Massot (Nov 22 2019 at 16:36):

The slightly confusing rule is: rw rewrites every occurrence but it has to be the exact same rewriting. For instance, if your goal contains ((a+b)+ (c+d)) + (a + b), rw add_comm will perform only one rewrite, although there are five additions, because it will see the outermost addition (second to last one in left to right order), specialize the add_comm lemma to add_comm ((a+b)+(c+d)) (a+b) and perform all rewrite using that specialized lemma. On the other hand rw add_comm a b will perform two rewriting, because that specific lemma matches two subterms in the goal.

Last updated: May 17 2021 at 23:14 UTC


text:
## On simultaneous approximations by irreducible fractions.(Russian)Zbl 0464.10040

### MSC:

 11K60 Diophantine approximation in probabilistic number theory


text:
SMS scnews item created by Florica Cirstea at Thu 13 Aug 2009 1954
Type: Seminar
Distribution: World
Expiry: 17 Aug 2009
Calendar1: 17 Aug 2009 1500-1600
CalLoc1: Carslaw 454
Auth: florica@p7192.pc (assumed)

# PDE Seminar: Brandolini -- Some recent estimates on linear Neumann eigenvalues and eigenfunctions

 Speaker: Barbara Brandolini from University of Naples

3pm-4pm, Monday, 17 August 2009, Carslaw 454

Title: Some recent estimates on linear Neumann eigenvalues and eigenfunctions

Abstract: We prove a reverse Holder inequality for the first nontrivial eigenfunction of
the Neumann Laplacian.  From this we deduce a lower bound for the first nontrivial
eigenvalue in terms of the isoperimetric constant relative to the domain.  We finally
show that our results are asymptotically sharp at least in dimension two.

Before the talk, we will take the Speaker to lunch at the restaurant Thai@home on Glebe
Point Road.  We will meet at 1pm on the outdoor Plaza on the 2nd floor of Carslaw.

Actions:
Calendar (ICS file) download, for import into your favourite calendar application
UNCLUTTER for printing
AUTHENTICATE to mark the scnews item as read
School members may try to .


text:
# Constraint programming using solver: how to state the time limit of search?

Could someone give additional explanations about the TimeLimit of CPLEX for CPO?

I am using CPLEX to solve a problem of assignment and scheduling in which each task needs to be assigned to a worker, by respecting a correct sequence. There is a cost of changing the worker and the cost of tasks execution changes from one worker to another.

When I do not state the time limit, it takes about 15 hours to return a result. When I stated the time limit = 10 seconds, it took about 10 seconds to find a correct result, correspondent to that found by the solver after turning during 15 hours.

The only difference between both results is the sequence of tasks executed by the same worker, but it does not imply on the objective function.

My questions are:

1. Why the solver takes 15 hours to obtain a result if it could take just few seconds to obtain an equivalent result?
2. Could it be related to the fact that there are many possibilities of sequencing the tasks without changing the objective function value?
3. How to know which is the most appropriate time limit to use?
• You mention both CPLEX and CPO (CP Optimizer). Which one are you using? – prubin Dec 5 '19 at 18:56
• I mean I use the docplex.cp module to call the IBM ILOG CPLEX Optimization Studio in Python – campioni Dec 6 '19 at 9:31
• @campioni, would you say, what is your problem scale? (I mean, how many binary variables are there in your model?). AFAIK, CP optimizer logic is different from CPLEX logic to solve the MILP problems. I tried a couple of large scales scheduling problems using CP and they were solved in a reasonable time. Are you sure to use the correct form of the CP definition in your model? – A.Omidi Dec 6 '19 at 10:55
• @A.Omidi I used both solvers (CPLEX and CPO) but my model is based in mathematical programming, I mean I did not use specified functions in CPLEX, such alldifferent or even interval variables to state precedence. In think it spend a long time to be solved through CPO because of the model symmetry as they said in other answers. – campioni Dec 9 '19 at 14:21

What Nikaza said, and also: problems of this kind often have a large number of solutions that differ only by trivial permutations/relabellings/reorderings. For example, if one worker is assigned 10 tasks that can be done in any order with no change to the objective, then there are 10! = 3628800 solutions that will have exactly the same OF. If another worker has the same situation, then you have around 1013 optimal solutions that are identical modulo task reorderings, and so on.

In my experience, this kind of symmetry sometimes leads to excessively slow solution times with the sort of behaviour you describe, and tweaking the problem to break symmetry can significantly improve performance.

A couple of options for symmetry-breaking:

• Add constraints that impose orderings.

For example, if tasks A and B are identical in requirements, we might require that B can't start before A does.

• Make small perturbations to the objective function weights such that the solutions are no longer identical.

For example, if your model has a binary decision variable x[j,t] indicating whether task j begins at time t, then adding the sum of x[j,t]*r[j,t]*smallnum to the objective function breaks symmetry with respect to task ordering.

Here, r[j,t] is a randomly-determined parameter chosen separately for each [j,t] and smallnum is chosen small enough that the solution will still be approximately optimal with respect to the original objective function, but large enough not to get rounded away to zero.)

IME, these options can sometimes produce a surprisingly large improvement in runtime.

• Thank you for your explanations. For your second suggestion, in practice, how could I make perturbations on the objective functions ? – campioni Dec 6 '19 at 16:33
• @campioni I added an example, though you may need to adapt to fit your model. – Geoffrey Brent Dec 6 '19 at 22:16
• thanks for your help, I appreciate it. – campioni Dec 9 '19 at 13:27

CPLEX is probably spending 15 hours trying to reduce the optimality gap. This is very common issue in MILP problems.

The way MILPs are solved is commonly to make all integers continuous and solve the resulting LP. That is called the linear relaxation, and its solution provides a valid lower bound (if we are minimising) on the objective function. The difference between the best integer solution the solver has found and that lower bound is called the optimality gap.

In branch and bound, it's common to find the best integer solution fairly quickly, and the bulk of calculations is spent to prove that this solution is globally optimal, i.e., to close the optimality gap.

There are two ways to handle that in practice:

1. Use the optimality gap as your convergence criterion, i.e., "stop calculating if the gap is smaller than 10%".
2. Define a design goal in your objective function. This means that, instead of minimising cost you would minimise cost-good_enough_cost and stop when that number is zero. This could also be expressed in terms of a constraint: cost-good_enough_cost >= 0.
• nikaza, thanks for replying. I am trying to optimize a 0-1 integer linear programming problem. In this case, optimality gap could also be an issue? – campioni Dec 5 '19 at 13:09
• Yep. This is how CPLEX is likely to solve it under the hood, it will relax the binaries to continuous between 0-1 and solve a series of LPs in branch-and-bound to get lower bounds (among other tricks). – Nikos Kazazakis Dec 5 '19 at 16:27
• Thanks for your explanations. In docplex.cp, the parameter to state this optimality gap is 'OptimalityTolerance' ? – campioni Dec 6 '19 at 9:46
• Just for the record, I think the poster is using CP Optimizer rather than CPLEX, so gap reduction is not exactly what is going on. – prubin Dec 7 '19 at 16:19
• @prubin, at first, thanks for your explanations, it has helped me a lot. I've solved with CPLEX and CPO (TimeLimit=10seconds) and I compared both results. Fixing the TimeLimit, CPO is much more faster than CPLEX. I tried to state the gap tolerance in CPLEX, but it did not help to increase the search speed significantly. – campioni Dec 9 '19 at 13:50

CP Optimizer is doing a search for progressively better solutions. Each time it finds a feasible solution, it restricts the remainder of the search to solutions with better objective values than what it just found. So you should get a sequence of progressively better solutions. Given enough time and memory, it will exhaust the search space, and the last solution found will be provably optimal. Stopped early, the last solution found will not necessarily be optimal.

So it sounds as though CPO is finding a good solution in the first 10 seconds and then struggling to find a better one. If the 15 hour run completed the search (i.e., was not stopped by time or memory limits), then that early solution is optimal.

As Geoffrey Brent said, you might have significant symmetry in your model, in which case symmetry-breaking constraints would reduce the size of the search space (and hopefully let CPO terminate sooner). There are also ways you can change the search strategy CPO uses, but if the early solution is actually optimal, I'm not sure that messing with the search strategy will get proven optimality (exhaustion of the search space) any sooner.


text:
• Subject:

...

• Topic:

...

Two wires of copper having the length in the ratio 4 : 1 and their radii ratio as 1 : 4 are stretched by the same force. The ratio of longitudinal strain in the two will be

(a) 1 : 16                                             (b) 16 : 1

(c) 1 : 64                                             (d) 64 : 1